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1.0  Summary 

The  overall  technical  goal  of  this  project  is  to  develop,  test  and  validate  a  novel  code  for  Large- 
Eddy  Simulation  (LES)  of  reacting  flows  at  high  speeds.  The  code  consists  of  two  main 
components: 

1 .  Flow  solver:  The  flow  solver  is  based  on  a  high-order,  Discrete  Spectral  Element  Method 
(DSEM)  on  unstructured  grids.  Shocks  and  other  discontinuities  are  simulated  using  a 
combination  of  Entropy  Viscosity  (EV)  method  and  spectral  filtering.  LES  subgrid  scale 
effects  are  accounted  for  implicitly  using  explicit  filtering  with  negligible  computational 
overhead. 

2.  Combustion:  Chemical  reaction  is  simulated  using  Filtered  Mass  Density  Function 
(FMDF)  method.  This  method  directly  simulates  the  species  transport  and  provides  a 
closed  fonn  for  the  reaction  source  term.  Thermophysical  properties  are  calculated  using 
in-house  routines  developed  under  this  project. 

Full  integration  of  FMDF  into  the  DSEM  code  is  planned  for  Phase  II;  however,  the  integration 
process  began  in  Phase  I.  During  Phase  I,  a  compact  Finite  Difference  (FD)  code  was  used  for 
testing  various  aspects  of  FMDF. 

This  Final  Report  is  composed  of  three  main  components.  The  first  component  is  the 
implementation  of  the  EV  method  in  the  DSEM  code  for  ID  scalar  transport  and  2D  Euler 
equations.  The  second  component  deals  with  development  of  FMDF  routines  for  supersonic 
flow,  and  includes  details  of  calculation  of  thermophysical  properties  and  the  results  for 
hydrogen-air  combustion  in  mixing  layers.  The  third  component  deals  with  implementation  of 
FMDF  in  the  DSEM  code. 

In  Phase  I,  proof  of  concept  was  provided  in  all  of  the  three  aforementioned  components.  In  the 
first  component,  the  EV  method  was  extended  to  DSEM  for  both  ID  scalar  transport  and  2D 
Euler  equations.  Several  tests  were  conducted  to  show  the  effectiveness  of  EV  in  capturing 
shocks  and  removing  spurious  oscillations.  Spectral  filtering,  however,  proved  not  to  be  effective 
in  DSEM.  In  the  second  component,  several  routines  were  developed  for  calculation  of 
thennophysical  properties  of  the  reactive  species  and  the  mixture.  These  routines  were 
successfully  validated  in  detailed  zero-dimensional  flow  for  the  stirred  reactor.  Then,  several 
simulations  were  conducted  for  combustion  of  hydrogen-air  in  two-  and  three-dimensional 
mixing  layers  for  both  subsonic  and  supersonic  flows  to  demonstrate  the  capabilities  of  the  new 
FMDF  routines.  For  the  third  component,  the  focus  in  Phase  I  was  on  implementation  of  particle 
tracking  algorithm  and  statistical  averaging  scheme  in  the  DSEM  code.  A  consistency  test  was 
conducted  for  the  fluid  density  in  a  free  mixing  layer  and  it  was  shown  that  the  Lagrangian 
(FMDF)  and  Eulerian  (DSEM)  produce  consistent  results  for  density. 
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2.0  Introduction 


In  the  past  few  years,  several  different  groups  have  invested  significant  efforts  to  develop 
hypersonic  vehicles.  In  May  of  2010,  the  news  media  reported  the  longest-ever  supersonic 
combustion  ramjet-powered  flight.  The  200+  second  burn  by  the  X-5 1  ’s  Pratt  &  Whitney 
Rocketdyne -built  air  breading  scramjet  engine  accelerated  the  vehicle  to  Mach  5.  The  program 
had  a  goal  for  a  300  s  flight  to  Mach  6  but  after  140  s  of  operation  some  anomalies  were  noticed, 
requiring  the  controllers  to  activate  the  self-destruct  function.  While  this  is  a  major  step  to  make 
supersonic  combustion  a  reality,  it  also  brings  to  focus  the  significant  challenges  that  remain 
ahead  in  designing  supersonic  combustion  engines.  These  tests  are  very  expensive  to  conduct 
and  require  a  long  period  of  preparation.  An  alternative  approach  could  be  provided  by  modeling 
and  simulation  but  the  existing  models  are  too  pre -mature  for  tackling  this  problem. 

Computational  Fluid  Dynamics  (CFD)  is  rapidly  expanding  as  a  viable  design  tool  in  many 
applications  in  order  to  minimize/avoid  costly  experiments.  Nevertheless,  the  rate  of  penetration 
of  CFD  into  the  design  sequence  has  not  been  the  same  in  all  areas.  One  area  of  particular 
interest,  and  yet  extremely  challenging,  involves  combustion  systems,  e.g.  internal  combustion 
engines,  gas  turbines,  and  combustors.  The  main  challenge  in  these  systems  arises  from  the 
presence  of  turbulence,  combustion,  and  (in  some  cases)  a  dispersed  phase  of  droplets  or  solid 
particulates.  Each  of  these  phenomena  is  characterized  by  the  existence  of  a  wide  range  of  scales, 
and  exhibits  tremendous  difficulty  from  a  modeling/simulation  standpoint.  Consequently, 
accurate  modeling  of  these  phenomena  has  remained  an  unresolved  issue  after  decades  of 
intensive  research. 

The  situation  is  even  more  complicated  in  supersonic  combustion  modeling  (e.g.  in  scramjets 
and  afterburners)  since  the  flow  residence  time  in  the  combustor  is  considerably  short,  in  the 
order  of  milliseconds.  The  high  speed  of  the  flow  significantly  affects  the  nature  of  turbulence- 
combustion  interactions  by  imposing  flow  time  scales  comparable  to  those  of  combustion. 
Traditional  models  built  upon  the  assumption  of  fast  chemistry  (as  compared  to  the  flow  time 
scale)  may  no  longer  be  valid  for  supersonic  combustion  and  new  models  are  needed  to  take  into 
account  the  comparable  time  scales  of  the  flow  and  combustion. 

When  the  reacting  flow  furthermore  contains  shocks,  the  scale  range  of  the  continuum  flow  is 
enonnous.  The  tremendous  complexity  of  the  problem  has  posed  the  highest  of  demands  on 
numerical  methods  and  models  and  has  left  many  physical  situations  not  fully  understood. 
Improved  understanding  of  shock-turbulence  interaction,  and  control  of  shocked  flows  is 
eminent  to  improving  combustor  efficiencies,  supersonic  devices,  and  high-speed  flow 
environments. 

The  overall  focus  of  this  proposed  project  is  on  modeling  and  simulation  of  high-speed  turbulent 
reacting  flows.  The  Enabling  Energy  Systems  (EES)  Inc.  has  assembled  a  team  of  highly 
experienced  experts  to  tackle  all  the  main  issues  involved  in  modeling  of  such  complex  flows. 
The  members  of  this  team  are  proposing  several  innovative  ideas  and  have  a  long  history  of 
working  together  with  complementary  expertise.  The  ultimate  goal  of  the  project  is  to  develop 
advanced  CFD  software  based  on  a  high-order  spectral  element  method  on  unstructured  grids. 
The  combustion  modeling  will  be  via  the  accurate,  pdf-based,  FMDF  method  that  we  have 
developed  for  the  LES  of  turbulent  reacting  flows.  Shock  discontinuities  will  be  captured  using 
an  entropy  viscosity  approach. 
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In  the  next  section,  we  explain  methods,  assumptions,  and  procedures  for  each  of  the  three 
components  of  this  project.  This  is  followed  by  presentation  and  discussion  of  the  results. 
Finally,  some  concluding  remarks  are  provided  in  the  last  section. 
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3.0  Methods,  Assumptions,  and  Procedures 

3.1  Implementation  of  the  EV  Method  for  1 D  Scalar  Equation 

In  this  section  we  consider  the  one-dimensional  linear  and  nonlinear  equations  for  a  scalar  in 
order  to  describe  the  implementation  of  the  EV  method  into  the  DSEM  code.  The  simplicity  of 
this  problem  allows  us  to  develop  a  better  understanding  of  the  method.  We  then  extend  this 
implementation  to  the  two-dimensional  Euler  equations  in  the  next  section. 

3.1.1  Mathematical  Modeling 

We  briefly  present  the  mathematical  background  of  the  DSEM  method  for  a  one-dimensional 
scalar  conservation  equation  as  a  relevant  basis  for  describing  explicit  filtering  and  entropy 
viscosity  methods.  For  a  detailed  description  of  the  DSEM  method,  we  refer  to  [1,2]. 

Governing  Equations 

We  consider  the  one-dimensional  scalar  conservation  equation  with  the  relevant  initial  and 
boundary  conditions, 


+/x(w)  =  0, 

x  G  [a,  b],  t  >  0 

(1) 

it(x,  0)  =  it0(x) 

(2) 

u{a,t )  =  gL(t), 

u{b,t)  =  gR(t ) 

(3) 

Here,  we  assume  the  flux  has  two  components:  advective  flux  and  diffusive  flux  that  comes  from 
the  EV  method,  f  —  fA  +  fEV  ■  In  this  project,  we  consider  the  linear  advection  and  the 
nonlinear  Burgers  problems.  These  two  problems  can  be  viewed  as  two  special  cases  of  the 
general  Equation  (1).  If  we  set  f  —  u  and  /  =  ^  u 2  in  Equation  (1)  we  obtain  the  linear 
advection  and  the  nonlinear  Burgers  equations,  respectively. 

Discontinuous  Spectral  Element  Method 

The  DSEM  is  a  collocation  method  for  the  solution  of  compressible  flows  on  staggered  grids. 

The  solution  values  are  defined  at  the  nodes  of  a  Gauss  quadrature  rule,  and  the  fluxes  are 
computed  at  the  nodes  of  a  Gauss-Lobatto  rule.  This  method  is  conservative,  a  feature  that  is 
desirable  for  application  of  shock  capturing  techniques.  Figure  1  shows  a  representation  of 
Gauss-Gauss  points  and  Gauss-Lobatto  points  in  a  ID  domain. 

■-0-0 - O - O - O - O - O - •— O — 0 

Figure  1:  Staggered  Arrangement  of  Solution  Variable  and  Fluxes 
Open  Circles:  Gauss-Gauss  Quadrature  Points,  Solid  Squares:  Gauss-Lobatto  Quadrature  Points 

In  this  part  of  the  project,  we  use  the  ID  DESM  code.  The  domain  [a,  b]  is  subdivided  into 
multiple  non-overlaping  subdomains  Ok  —  [ak,  bk],  k  —  0,1, ... ,  K,  which  are  ordered  left  to 
right.  A  linear  mapping  is  used  to  map  the  subintervals  to  a  unit  interval.  Under  this 
transformation,  Equation  (1)  can  be  written  on  each  subinterval  as 
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Ut+fx(u)  =  0 


(4) 


Where  u  —  xxu  and  xxis  the  mapping  operator.  Below,  we  describe  our  numerical  method  based 
on  the  mapped  Equation  (4);  however,  we  drop  the  A  for  brevity. 

The  computational  domain  (Q)  is  represented  by  the  union  of  non-overlapping  elements,  Dk, 

n  =  (5) 

In  each  subdomain,  the  values  of  the  solution  or  fluxes  can  be  defined  on  Gauss  and  Gauss- 
Lobatto  points: 

Q(x)  =  ZUQMx) 

N-l 

Qw = y  o  iii  i® 

/  *  J  '  9  J  '9 
7=0 


(6) 

(V) 


These  polynomials  interpolate  values  of  any  function  Qj  on  the  Lobatto  grid  and  the  value  of 
Qj+rn  on  the  Gauss  grid.  lj  and  hj+m  are  the  Lagrange  interpolating  polynomials  of  order  N  and 
N-l  defined  on  Lobatto  and  Gauss  grid  points,  respectively, 


fy+1/2 


S-Xi 


1+1/2 


Xj  +  1/2  Xi+ 1/2 


i*j 


(8) 

(9) 


The  Gauss  quadrature  points  and  Gauss-Lobatto  quadrature  points  are,  respectively,  defined  as 


Element  Level  Filtering 

In  addition  to  the  EV  method,  we  also  consider  explicit  filtering  as  a  part  of  our  approach  to 
capture  shocks  and  stabilize  the  solution  while  preserving  high-order  resolution.  We  employ  an 
element  based  filtering  approach  following  Blackburn  and  Schmidt  [3].  Since  our  polynomial 
basis  functions  do  not  form  a  hierarchical  set,  filtering  can  be  implemented  by  projecting  to  a 
lower-order  set  of  basis  functions  in  the  same  family.  Defining  ijf  as  the  operator  that 
interpolates  a  polynomial  of  order  N  with  Np  =  N  +  1  nodes  onto  a  set  of  Mp  =  M  +  1  nodal 
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points,  we  produce  the  Lagrange-interpolant  projector  F  =  ijjljj1.  This  projects  the  original  nodal 
values  Q(xk)  to  another  set  at  the  same  nodal  locations  xk,  but  derived  from  the  Lagrange 
interpolants  through  the  smaller  (M+  1)  set  of  Gauss-Lobatto-Lobatto  nodes.  The  operators 
In  and  1$ can  be  derived  in  any  appropriate  fashion  from  the  unique  polynomials  passing  through 
the  corresponding  points. 

The  solution  of  the  discretized  equations  is  collocated  on  nodal  points  in  physical  space.  If  lj  are 
the  polynomial  basis  functions  then  the  nodal  values  of  any  smooth  function  Q(x)  could  be 
represented  by  the  spectral  expansion, 

N 

Q(*)  =  ^Q;W  (12) 

7=0 

By  applying  the  filtering  described  above  we  first  project  the  solution  to  a  lower  polynomial 
space  and  then  interpolate  them  back  to  higher  polynomial  space, 

Q,1lt=/w  QIm  (13) 

Entropy  Viscosity  Method 

The  entropy  viscosity  method  used  in  this  work  follows  the  work  of  Guermond  et  al.  [4] 
presenting  a  new  method  to  avoid  spurious  numerical  oscillations  in  high-order  numerical 
approximations  for  nonlinear  conservation  laws.  This  method  is  founded  on  the  idea  of  adding  an 
artificial  nonlinear  dissipation  tenn  to  the  numerical  discretization  as  an  alternative  to  previous 
methods  like  using  limiters  or  non-oscillatory  reconstruction.  The  magnitude  of  the  dissipative 
tenn  is  detennined  based  on  the  local  entropy  production  rate.  Since  large  entropy  production 
takes  place  near  shocks,  the  entropy  viscosity  tenn  is  large  in  this  region  and  small  elsewhere. 

In  this  report,  we  focus  on  scalar  conservation  equations,  which  have  many  entropy  pairs  and  we 
can  find  at  least  one  entropy  function  that  satisfies  an  entropy  inequality.  The  entropy  satisfies  a 
conservation  equation  in  the  regions  where  the  solution  is  not  oscillatory,  and  where  we 
experience  oscillations  or  shocks  the  entropy  production  will  grow  and  so  does  artificial 
dissipation  term  which  helps  to  damp  the  oscillations  and  obtain  a  smooth  solution  in  these 
regions  of  the  flow.  For  a  scalar  conservation  equation  we  assume  f2  to  be  an  open  connected 
domain  and  consider  a  model  problem  in  this  domain  as 

dtu(x,  t)  +  V.f(u(x,  t))  =  0,  x  £  at  >  0  (14) 

This  problem  is  solved  subject  to  an  initial  condition  and  proper  boundary  conditions.  It  is  well 
known  that  the  initial  boundary  value  problem  has  a  unique  entropy  solution,  which  satisfies  the 
differential  equation 

dtE(u)  +  V. F(u)  <  0,  xE  at  >  0  (15) 

for  any  pairs  E  and  F  such  that  E  is  convex  and  F(u )  =  /  E'(u)  fr(u)du.  The  function  E  is 
called  entropy  and  F  is  the  associated  entropy  flux.  It  is  proven  that  for  one  space  dimension 
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E(u )  =  —  yields  a  unique  entropy  solution.  If  we  consider  an  entropy  pair  (E,  F),  we  define  the 
entropy  residual  as 

Dh(x,  t)  =  dtE(uh(x,  t))  +  V.  F(uh(x,  t)),  xE  EL,t>  0  (16) 

And  then  we  define  a  viscosity,  vE,  based  on  this  residual 

vE(x,t)  =  cEh2(x)R(Dh(x,t))/\\E(uh(x,t))  -  E(uh(x,t))\\a  (IV) 

where  h(x)  is  the  local  mesh  size  at  x;  E  (iih  (x,  t))  is  the  space  averaged  value  of  entropy;  cE  is 
a  tunable  coefficient  which  depends  on  the  spatial  discretization  and  also  the  specific  problem  at 
hand.  The  simplest  function  one  can  consider  for  R  is  R(Dh )  =  \Dh  \  to  prevent  negative  values 
for  the  residual.  We  also  introduce  an  upper  limit  for  the  artificial  viscosity 

Vmaxfa’  0  ^max^-max  (A)  \f  (u(x,  t))  (18) 

Here  T  is  a  neighborhood  of  x  and  f'(u(x,  t))  is  the  local  wave  speed.  Then  we  define  the 
artificial  entropy  viscosity  as  follows: 

vh  =  S(minCiW,v£))  (19) 

where  S  is  a  smoothing  operator.  The  whole  procedure  is  equivalent  of  adding  an  artificial 
dissipation  term  —V.  ( vhVuh )  to  the  discritized  equations. 

3.1.2  Numerical  Methodology 

This  section  briefly  describes  the  implementation  of  element  level  filtering  and  entropy  viscosity 
method  in  the  DSEM  code. 

Element  Level  Filtering 

An  explicit,  low-pass  filter  is  used  for  filtering  operation.  Spectral  filtering  can  be  constructed 
using  either  discrete  polynomial  transform  (DPT)  or  interpolant-projection  (see  [3])  over  each 
element.  DPT  filtering  can  be  conveniently  applied  for  methods  with  modal  basis.  For  methods 
with  nodal  basis,  the  solution  has  to  be  first  transformed  to  modal  basis  before  the  DPT  filter  can 
be  applied.  Projection  filtering  on  the  other  hand  can  be  constructed  directly  on  the  nodal  basis. 
Since  it  does  not  require  an  extra  transformation,  interpolant-projection  filtering  is  more  efficient 
than  DPT  for  methods  with  nodal  basis.  Therefore,  for  our  nodal  basis  we  use  an  interpolant- 
projection  filter. 

In  the  interpolant-projection  filtering  procedure,  the  filtered  variable  of  degree  N  is  obtained  by 
projecting  the  variable  back  and  forth  to  a  lower  order  approximation  of  degree  M  defined  on  a 
subset  of  the  original  nodal  values.  As  a  first  step,  the  original  function  is  interpolated  from  a 
polynomial  degree  N  to  a  polynomial  of  lower  degree  M, 

N 

Q'(xi)=^lj(xi)Q(xj)  (20) 

j=o 
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where  xi?  Xj  are  the  nodes  corresponding  to  Pm  and  Pn,  respectively,  lj  G  Pn  is  the  Lagrange 
interpolating  polynomial.  The  above  operation  can  be  cast  in  terms  of  matrix-vector  product, 


QrTQj 

where. 


(21) 


N 


1?=  n  i=o’---’M’  j=°  - 

A,-lvxrxk 


,N 


k=0,k^j 


In  the  second  step,  the  function  Q '  (x)  is  projected  back  to  the  polynomial  space  N  giving  the 
fdtered  function, 


N 


Qfllt(xe)=^  lf<xe)Q'(Xf) 
^0 


(22) 


Again,  the  above  can  be  cast  in  matrix-vector  form, 

nnu  ,pr°n/ 
ve  Aef  Vf 


(23) 


where 


Tpr°  _ 

Aef 


V 


n 

k=0,k# 


xe  xk 
xf-xk  ’ 


e=0,....,N, 


f=0...,M 


In  the  staggered  grid  method,  this  interpolation-projection  operation  could  be  applied  to  both  the 
nodal  sets  (Gauss-Gauss  and  Gauss-Lobatto  nodes).  We  apply  the  filter  on  the  Gauss-Lobatto 
basis  since  it  preserves  the  end  values  of  the  original  function  and  ensures  C°  continuity. 

Entropy  Viscosity  Method 


In  the  spectral  element  method  in  one  space  dimension,  the  domain  is  divided  into  non¬ 
overlapping  segments.  Each  subdomain  K  is  an  image  by  a  map  gicof  the  reference  element  K 
=(-1,1).  We  define  the  quantity  h(x)  to  be  the  distance  between  two  consecutive  collocation 
points.  For  one  space  dimension,  the  defined  parameters  in  the  previous  section  will  take  the 
following  forms: 


h(xi)  =  min  xt  —  Xj  i  A  j 


(24) 


R(.Dh(Xi,  t))  =  |Z4(*i,t)| 


(25) 


vE(Xi,t)  =  cEh2(Xi ) 


_ I  DhXXy  01 _ 

-E(uh(xut))\\a 


(26) 


In  the  above  equations,  X;  represents  an  arbitrary  point  in  one  subdomain  K.  The  maximum 
viscosity  is  defined  within  an  appropriate  neighborhood  Vx  —  K  3  x.  The  maximum  viscosity  is 
defined  on  each  element  as 
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Vmax  I K  =  CmaxhKmax\f'(.uh(xi,  t))|  ,XiEK 


(27) 


Our  numerical  experiments  show  that  without  using  smoothing  operator,  the  results  would  be 
oscillatory  especially  for  higher  polynomial  orders.  Smoothing  is  done  for  each  element  and  we 
use  the  following  smoothing  operator 

S(<p)  =^(<K*t-i)  +  <p(.xi+1)  +  2<p(*f))  (28) 

The  time  stepping  method  used  in  this  work  is  the  low-storage,  third-order,  Runge-Kutta  scheme 
presented  by  Williamson  [6].  The  implementation  requires  that  we  store  the  value  of  the  solution 
at  Gauss  points  for  two  time  steps  in  order  to  be  able  to  evaluate  the  entropy  residual, 

°h  =  2 ~  4^“ 1  +  E^2)  +  VE*  (29) 

In  the  above  equation,  At  is  the  time  step  size  which  is  assumed  to  be  constant,  E £  = 

E(u %)  and  E^1  —  E(u]!L~1')  .  The  residual  is  evaluated  explicitly  using  second-order  backward 
difference  for  evaluating  time  derivative.  Having  evaluated  the  entropy  residual,  we  can  evaluate 
viscosities  vE  and  vmax  and  eventually  vh  which  all  are  dependent  also  on  the  tunable 
coefficients  CE  and  Cmax.  The  values  of  CE  and  Cmax  are  tuned  by  starting  from  a  large  value  for 
CE  so  that  only  the  first  order  viscosity  is  active,  i.e.  vh  —  vmax.  Then  we  reduce  the  value  of 
Cmax  until  the  method  is  on  the  edge  of  loosing  stability  or  smoothness. 

Algorithm  in  One  Space  Dimension 

Summarizing  what  have  been  described  so  far  about  staggered-grid,  multidomain  spectral 
method,  element-level  filtering  and  entropy  viscosity  methods  and  combining  them  we  arrive  at 
the  following  algorithm  for  approximating  the  scalar  problem  in  Equation  (1)  in  one  space 
dimension  with  explicit  filtering  and  entropy  viscosity: 

Algorithm  for  Staggered-grid,  Scalar,  ID 

1.  Use  solution  at  Gauss  points  to  construct  solutions  at  the  Lobatto  points 

u«  =  =  hft  1/zOi) 

2.  Compute  the  advective  flux: 

a.  Compute  the  flux  values  at  the  internal  point  at  Lobatto  grid 

F]A'K=fa(Uf ) 

b.  Apply  boundary  conditions 

F0A'°=fa(gL) 

Fo’K=fai9R) 

c.  Apply  interface  conditions 

•  If  the  wave  is  right  running  characteristics 

pA,K-l_pA,K 

Gv  -ro 

•  If  the  wave  is  left  running  characteristics 

pA,K  _  pA,K-l 
^0  —  tN 

3.  Compute  entropy  viscosity  diffusive  flux: 

a.  Average  solution  at  the  interfaces 
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u^u*-1  =  \{Uk  +  ut1) 

b.  Construct  entropy  on  Lobatto  points 

E?  =  \(Uft2 

c.  Calculate  average  entropy  over  the  domain 


h.+ 1 


K 


(ES(x,))  +  Ei;(x,)))/^htH 


l  l 

d.  Calculate  maximum  local  entropy 

maxE  —  max  (ijf),  i  =  l....N,  k  —  1 . K 

e.  Construct  entropy  time  derivative  on  Lobbato  points 

Af  1 00  =  2^0 EZ(x)  -  4 EZ~\x)  +  ^(x)) 

f.  Construct  entropy  residual  on  Lobatto  points 

/I 


Dfrixi)  —  max 


A  E(xf)  / 

'  (<Oi+i)  -/'(<(Xi))) 

2A  t 

h.  i 
t+2 

V 


A^(xj+1)  / 

+ 

H- 

1 

/?s' 

k 

2At 

h-i+i/2 

/ 

where,  (h.+i)  =  |*i+i  ~  xt\K 
g.  Calculate  entropy  viscosity 

Vhfri )  =  max  f  Cmaxh  i,  CEh2  1D^  (xj/  max  E 


t+2 


t+2 


h.  Differentiate  and  compute  diffusive  flux  at  Gauss  points 

pEV.K  _  y^Df/fe(D);j  =  l'j(Xi  +  1/2) 

i.  Average  flux  values  at  the  interfaces 


=F, 


w 


+  E, 


N 


L) 


4.  Combine  fluxes  to  get  the  total  flux 

Ff  =  F/’fc  +  FjEV,kj  =  0,1 . N-  k  =  1,2, ... .  K 

5.  Differentiate  flux  and  evaluate  the  solution  on  Gauss  grid  using  RK3 

jt  Uk+ 1/2  +  (DFk)j+1/ 2=0  j=0,l N-l  ;  k  =  1,2 . K 

6.  Filter  the  solution 

a.  Construct  interpolation  and  projection  matrices 


ljjlt-  FT  i=0 . N,  j=0..., 

1  1  I  Xj-Xi, 

k=0,k*j 
M 

■n 


Nfiit 


Tpro 

!ef 


k=0,k*f 


XfXk 


f=0 . N,  e=0 . N. 


Jilt’ 
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b.  Interpolate  solution  from  Gauss  to  Lobatto  Points 

=  /£t7K(/£)u  =  fi/+1/2Oi) 

Interpolate  from  Lobbato  points  (N)  to  lower  polynomial  space  (Nmt) 


c. 


uL= iintuk 


d.  Project  back  to  higher  polynomial  space 


U 


proj- 


Tpronfc 
1  uint 


e. 


Interpolate  from  Lobbato  to  Gauss  points 
Vftt,  =  lLUprolO!!)l.i  =  if  (*i) 

3.2  Implementation  of  the  EV  Method  for  2D  Euler  Equations 

As  the  second  step,  we  have  implemented  the  EV  method  in  the  two-dimensional,  inviscid  Euler 
solver.  Specific  accomplishments  include: 

1.  Development  and  testing  of  element  level,  explicit  low-pass  filtering  routines  in  2D 
DSEM  code  for  Euler  equations. 

2.  Development  of  a  consistent  and  compatible  numerical  formulation  of  the  EV  method  for 
2D  DSEM. 

3.  Programming  and  testing  of  2D  EV  subroutines  and  modules  and  integration  of  the  2D 
routines  into  the  2D  DSEM  inviscid  Euler  solver. 

4.  Preliminary  testing  of  the  2D  code  against  simple  benchmarks  for  supersonic  flows  with 
shocks. 


3.2.1  Mathematical  Modeling 

To  set  the  stage  for  the  discussion  on  the  implementation  of  the  EV  method  for  the  DSEM  based 
Euler  solver,  we  first  provide  the  conservation  laws  in  2D  and  discuss  filtering  methods  that  may 
be  used  in  combination  with  the  EV  method. 

Governing  Equations 

Since  the  EV  method  is  based  on  viscous  terms  that  have  the  same  form  as  the  physical  viscous 
stress  terms,  we  consider  the  governing  Navier-Stokes  equations  (with  viscous  fluxes).  The  non- 
dimensional  equations  for  conservation  of  mass,  momentum  and  energy  in  2D  Cartesian 
coordinates  are  given  by 


Qt  +  F 


Gy  = 


1 

R.6f 


(f;  +  op 


(30) 


where, 


Q  = 


Pi 

pu 

Fa  = 

>  1  X 

p 

p  +  pu2 

,  Gy  = 

pv  - 
puv 

pv 

puv 

p  +  pv2 

-pe. 

u(p  -1-  pe) 

v(p  +  pe). 

The  viscous  fluxes  take  the  following  form 
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(31) 


0 

0 

T11 

X21 

X12 

,  GX  = 

X22 

1 

'  y 

1 

pTu  +  VX12  +  (y  _  l)M2prTxJ 

|UT21+VT22  +  ^_i)Mf2prTyj 

where 

Tn  =  2[ux-  (ux  +  vy)/2],  X22  =  2[vy  -  (ux  +  Vy)/2],  Tn  =  x21  =  p(vx  +  uy) 

are  the  viscous  stresses  with  p  and  k  denoting  non-dimensional  temperature-dependent  viscosity 
and  conductivity  coefficients,  respectively.  In  the  EV  method,  the  physical  viscosity,  p,  is 
replaced  by  an  artificial  viscosity  that  damps  the  numerical  Gibbs  oscillations  induced  by  the 
high-order  polynomial  approximation  near  the  discontinuities. 

The  energy  flux  is  given  by 

Pe  =  +  P(u2  +  v2)/2  (32) 


In  the  non-dimensional  form  of  the  Navier-Stokes  equations  the  non-dimensional  Reynolds 
number,  Prandtl  number  and  Mach  number  appear  as 

Ref  =  U*fL}p}/p*;  Pr  =  C-^f-,  Mf  = 

where  *  denotes  the  dimensional  variables.  To  close  the  set  of  governing  equations,  the  ideal  gas 
equation  of  state  in  non-dimensional  form  is  given  by 


P  = 


PT 

yMf 


(33) 


Element  Level  Filtering 

Element  level  filtering  can  be  used  to  stabilize  spectral  simulations  of  shocked  flows.  Filtering 
may  be  employed  in  combination  with  EV.  We  are  considering  the  weak  filter  developed  by 
Blackburn  and  Schmidt  [3],  since  it  has  an  excellent  compatibility  with  our  nodal  DSEM 
formulation.  This  filter  projects  the  higher-order  solution  approximation  to  a  lower-order 
interpolant,  and  removes  the  high  frequency  content  of  the  high-order  interpolant.  Consequently, 
the  lower-order  intepolant  is  projected  back  to  the  original  nodes  of  the  higher-order  interpolant. 
Introducing  ijf  as  the  operatorthat  interpolates  a  polynomial  of  order  N  with  Np  =  N  +  1  nodes 
onto  a  set  of  Mp  =  M  +  1  nodal  points,  the  projection  is  summarized  as  F  =  1^1”,  with 


r M  _  {II?=o(*7  Xp)} / 

Njk  / ni=0(xk-xq ) 


j  =  0,  ...,M  k  =  0, ...  ,N 


(34) 
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where  the  x  values  are  the  coordinate  for  Gauss-Lobatto-Lobatto  (gll)  nodes.  Figure  2  shows  a 
schematic  of  the  two-stage  projection. 


2  r- 


Project 


1  b  ‘  I- 

0  h  oh 


I  1  I  1  I  i  I  I  1  i  I  1  1  i  1 


LLLLLLLUJ 


-l 


0 

x 


-1 


0 

x 


Figure  2:  Element-Level  Filtering  Through  a  Two-stage  Projection 

Courtesy  of  |3] 


Since  the  grids  in  2D  are  orthogonal  the  filtering  may  be  perfonned  along  ID  lines  in  both 
directions.  Effectively,  the  filtering  is  reduced  to  a  matrix-vector  multiplication  as  follows 

u=FuFt  (35) 


Where  F  is  a  one-dimensional  filter/projection  matrix. 

Entropy  Viscosity  Method  for  Euler  Equations  in  2D 

In  the  EV  method  an  artificial  entropy  viscosity  replaces  the  physical  viscosity  and  the  thermal 
conductivity  as 


kh  = 


Pr 

7-1 


Phi 


Ph 

vh=  — 

Ph 


The  entropy  viscosity  has  the  same  dissipative  effect  on  flows  as  the  viscous  stresses  in  the 
Navier-Stokes  equations.  The  entropy  viscosity  dissipates  the  unphysical  numerical  Gibbs 
oscillations  (energy)  in  higher-order  approximations  of  the  solution  generated  by  shock 
discontinuities.  While  it  is  desired  to  damp  the  numerical  oscillations,  physical  small-scale 
features  should  not  be  affected  by  the  artificial  damping.  To  ensure  that  only  high-frequency, 
numerical  oscillations  are  dissipated  by  the  entropy  viscosity,  but  not  the  small-scale  flow 
features  in  smooth  flow  regions,  the  entropy  viscosity  is  scaled  by  the  entropy  generation  in  the 
flow.  Shocks  generate  significant  entropy  and  hence  the  entropy  viscosity  is  large  in  the  vicinity 
of  shocks,  while  in  smooth  flows  the  entropy  generation  and  thus  the  entropy  viscosity  are 
significantly  lower. 


The  entropy  in  an  ideal  gas  flow  is  determined  as  [4] 


(36) 


The  entropy  generation  and  entropy  residual  are  detennined  as  follows: 

Dh(x,  t)  =  —  Sh(x,  t)  +  V.  ( uh{x ,  t)Sh(x,  t))  (37) 
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R{Dh(x,t ))  =  \Dh(x,t)\ 

The  entropy  viscosity  is  scaled  with  the  entropy  residual  as  follows: 
Pe O,  t)  =  CEph ( x ,  t) h(x)2  | Dh (x,  t)  | 


(38) 


(39) 


where  fiE  is  the  dynamic  viscosity.  In  the  EV  method,  we  define  another  similar  entropy 
generation  term  yielding, 


Dh,i  —  +  V.  ( uhSh ) 


Dh,2  Ph  Sh  (  Qf-Ph  hPh ) 


(43) 

(40) 


where  we  define  the  value  of  entropy  residual  as 
Dh  =  max(ph  l,  Dh  2) 

The  second  entropy  residual  Dh  2  essentially  accounts  for  inaccuracies  in  mass  conservation. 
Again,  similarly  to  ID  entropy  viscosity  method,  the  maximum  dynamic  viscosity,  fJ.max,  is 
evaluated  as 


M max  I K  Cmax  WPhWkIikWM  +yfyTh\\K 


(41) 


Finally,  we  set 

lih  =  min(iAmax,  fiE)  (42) 


3.2.2  Numerical  Methodology 

The  above  fonnulation  has  been  implemented  in  our  collocation  DSEM  approach  and  is 
summarized  as  follows: 


Algorithm  for  Staggered-grid,  Euler  Equations  in  2D 


1. 

2. 


At  Gauss  points  the  entropy  function,  S,  is  determined  based  on  the  functional 
relation  in  Equation  (36). 

The  entropy  residual  Dh  l  is  determined  as  follows: 

a.  The  physical  values  of  the  solution  from  Gauss-Gauss  (gg)  points  are 
interpolated  to  Lobatto-Gauss  (lg)  points  along  the  lines  in  both  grid 
directions  using  projection  matrices 

N- 1 N- 1 


m= 0  n= 0 


0 99 

V  1  1 

m+2’n+2 


N-l 


Qle(xuY  i)=  V  Q9\  xh  iTO 

V  1  +  2 '  vn-\-7j,TL-¥—  m  +  2 


m= 0 
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h  i  is  the  Lagrange  interpolating  polynomial  on  the  Gauss  grid. 

n  +  2 

The  values  of  the  interpolation  matrix  are  calculated  in  a  pre-processing 
step  and  stored  in  I. 

Similarly  in  y  direction 

In 


N-l 


Q9l{Xi+1/2,Yj)  =  + 

n= 0 


99 

,1  ,1 

m+2,n+2 


b. 


c. 


d. 


3. 


4. 


The  entropy  flux  values  ( uhSh )  are  computed  at  lg  points  from  the 
interpolants,  Q,  and  the  functional  relations 
Fx=uSh 
FyS=vSh 

A  Patching  process  is  applied  to  the  entropy  fluxes  so  that  we  make  sure 
the  entropy  fluxes  maintain  the  C°  continuity  over  boundaries  and 
interfaces. 

The  entropy  fluxes  are  differentiated  in  X  and  Y  directions  to  form  the 
divergence  of  the  entropy  flux  V.  ( uhSh ) 

%  =  2S  f“>  (x„  rj+ ,|)  DO',  OF1*  (xb  y/+|) 

%  =  Im  F“  (X1+1/2,  Yt)  aij(^1,2)=i:f.-o1  D(i.j)F,a  (xt,  F.+.) 

is  the  Lagrangian  interpolant.  The  derivative  comes  down  to  a  matrix- 
vector  product  D  *  F  with  D  a  differentiation  projection  matrix, 
e.  The  time  derivative  is  discretized  using  a  third  order  backward 
differencing  method. 

The  second  entropy  residual  Dh  2  is  determined  in  a  same  manner  as  Dhl  and 
find  Dh 

The  entropy  residual  from  Gauss  points  is  interpolated  to  Gauss-Lobatto  points 
to  determine  the  viscous  flux  terms. 


Interface  and  Boundary  Treatment 

Interpolation  of  the  solution  leads  to  different  solution  values  at  the  sub-domain  interface  points, 
one  from  each  of  the  contributing  sub-domains.  By  enforcing  continuity  of  the  advective, 
entropy  and  viscous  fluxes  over  the  interface,  the  elements  are  connected.  Moreover,  flux 
continuity  yields  a  conservative  method.  The  inviscid  fluxes  are  computed  using  an  approximate 
Riemann  solver.  Given  the  two  solution  states  1  and  Qq,  the  flux  in  each  spatial  direction  can 
be  expressed  as  [7]: 


Ca(Qw_1,Qo)  =  ^“(QJT1)  +  Fa(Qo)]  -  ^RAR-1^  - 


N 


A 


(43) 


where  Fa  is  the  vector  of  advective  fluxes  and  R  is  the  matrix  of  the  right  eigenvectors  of  the 
Jacobian  of  Fa  computed  using  Roe-average  of  Q/v  1  and  Qq.  At  the  boundaries,  the  physical 
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boundary  can  be  viewed  as  an  interface  between  the  external  state  and  the  computational  domain 
and  the  Riemann  solver  is  applied  between  external  specified  and  the  internal  solution  values.  A 
similar  approach  is  followed  for  entropy  fluxes.  We  used  a  Lax-Friedrichs  flux 

re(Qw-1,Qo)  =  ^[Fe(Q r1)  +  Fe(Qo)]  -  ^a(s(Qq)  -  S(Q^1))  (44) 

to  establish  continuity  of  the  entropy  fluxes  over  interfaces,  since  it  decouples  the  entropy 
approximation  entirely  from  the  flow  solver.  Moreover,  it  is  more  dissipative  than  the  Roe 
solver,  which  improves  the  stability  of  the  EV  method. 

3.3  FMDF  Model  Development  for  LES  of  High-Speed  Turbulent  Flows 

This  part  of  the  report  describes  our  efforts  on  the  development  and  testing  of  hybrid  LES/FMDF 
methodology  for  numerical  simulation  of  hydrogen-air  combustion  in  high  speed  turbulent  flows. 
As  we  mentioned  earlier,  in  Phase  I  a  compact  FD  code  was  used  to  test  the  FMDF  method.  To 
obtain  the  turbulent  gas  dynamic  field,  the  filtered  compressible  LES  equations  are  solved  with 
high  order  Eulerian  FD  methods.  The  hydrogen-air  mixing  and  combustion  are  implemented  by 
solving  the  FMDF  equations  via  stochastic  differential  equations  and  a  Lagrangian  Monte  Carlo 
(MC)  particle  method.  In  the  first  part  of  this  work,  different  sets  of  hydrogen-air  reaction 
models,  including  a  4-step  reduced  mechanism,  two  14-  and  16-step  skeletal  mechanisms,  and  a 
37-step  detailed  mechanism  are  successfully  implemented.  Different  methods  for  calculating  the 
thennophysical  properties  of  the  individual  species  and  the  mixture  are  also  investigated.  The 
most  efficient  and  accurate  methods  are  chosen  for  the  multi-species  calculations  of  scalar  and 
heat  transport  and  reaction  in  the  LES/FMDF  code  with  various  kinetics  models.  The  results  of 
the  reacting  and  non-reacting  compressible  turbulent  mixing  layers  are  discussed  in  Section  4. 
The  hybrid  LES-FD/FMDF-MC  method  is  employed  for  the  simulations  of  subsonic  and 
supersonic  reacting  and  non-reacting  mixing  layers.  In  this  report,  we  first  present  basics  of 
hydrogen-reaction  models  to  be  used  in  the  LES/FMDF.  We  will  then  describe  the 
implementation  of  these  reaction  models  into  LES/FMDF  code  and  the  results  obtained  by  this 
code  for  nonreacting  and  reacting  subsonic  and  supersonic  mixing  layers  are  presented. 

In  the  hybrid  LES/FMDF  methodology,  two  sets  of  Eulerian  and  Lagrangian  equations  are 
solved  together  for  velocity,  pressure  and  scalar  (temperature  and  mass  fraction)  fields.  These 
equations  are  presented  and  discussed  in  the  following  subsections. 

3.3.1  Filtered  Navier-Stokes  Equations 

The  filtered  Navier-Stokes  equations  including  continuity,  momentum  and  energy  are  expressed 
as  follows: 


dp  c{pu.) 

—  + - —  =  0 

dt  dx . 

J 


(45) 


d(pu.)  |  d(puu. )  _  gp  ^  da..  ^  dr., 

dt  dx .  dx.  dx .  dx . 

J  1  J  J 


(46) 
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(47) 


d(pet)  |  djpeu.) 
dt  dx. 

J 


d(  piij )  d(u.a..)  dq.  d[yRQ.I{y- 1)) 

dx.  dx.  dx .  dx . 

J  J  J  J 


Y(h°Wcb  ) 

v  a  a  a ' 


a= 1 


where  total  energy  et,  filtered  viscous  stress  tensor  and  heat  conduction  vector  are  defined  as, 


rr°  a= 1 


p  1 

—+-ukuk 

P  2 


(48) 


dT 

q  =  A(T,Y  )  — 

'  a’dx. 


(49) 


The  sub-grid  scale  (SGS)  stress  tensor  and  heat  flux  vector  z  and  Q  are  defined  respectively  as, 
t..  =  -p(uiuj  - u.u.)  (50) 


Q,=  P&u.-fuj) 


(51) 


In  order  to  close  the  governing  equations,  Smagorinsky  and  gradient-diffusion  models  are 
implemented.  Therefore,  the  modeled  SGS  stress  tensor  and  heat  flux,  r  and  Q,  respectively,  are 
written  in  the  following  forms: 


r.  - 

V 


=  Pot(Sg- 


(52) 


Pr  dx. 

‘  j 


where  eddy  viscosity  is  defined  as, 
ut=CsA2\S\  ,  ISH2SJ/2 


(53) 


(54) 


3.3.2  Compressible  Scalar  FMDF  Equations 

The  scalar  FMDF  represents  the  joint  PDF  of  the  scalar  vector  at  the  subgrid-level  and  is  defined 
as 


17 

Approved  for  public  release;  distribution  unlimited. 


(55) 


—  /•  +oo  —  — 

P.(x¥;x,t)=  [  p{x\t)3[^  ,0>{x\t)]G{x\x)  dx.' 

^  J  —oo 


(56) 


where  G  denotes  the  filter  function,  +  is  the  scalar  vector  in  the  sample  space,  and  <9  is  the 
“fine-grained”  density  [15].  The  scalar  vector  <t>  =  (pa,(a  =  1>— > Ns  +1)  includes  the  species  mass 

fractions  and  the  specific  enthalpy.  The  scalar  FMDF  transport  equation  is  obtained  from  the 
transport  equation  for  the  unfiltered  scalar  [16]: 


dpcpa  |  dpu.<pa 
dt  dxi 


_d 

dx. 


i  \ 


~S<P» 

dx, 


+p(s«+sr) 


i  J 


(57) 


Here,  for  simplicity,  we  consider  the  scalar  equation  in  the  Cartesian  coordinate  system.  For  the 
species  mass  fraction {ct  =  1  ,...,NS)  ,  the  source/sink  tenn  SR  =  coa  in  Eq.  (57)  represents  the 
production  or  consumption  of  species  a  due  to  the  chemical  reaction.  For  the  energy  or  enthalpy 

TV, 

0 a  =  (Vv+] )  ,  the  source  tenn  SR  =  -^  ( h°Wa(ba  (represents  the  heat  of  combustion,  and  the  tenn 

a= 1 


s 


cmp 

a 


1  dp 
—  (—  +  M, 
P  dt 


dp  du- 

—  +  <7 . . - 

dxt  dxj 


) 


is  due  to  compressibility  and  viscous  energy  dissipation.  The 


FMDF  transport  equation  is  obtained  by  inserting  the  instantaneous  unfiltered  scalar  equation 

&  (D 

into  the  time  derivative  of  fine-grained  density  —  = - — -  and,  filtering  that, 

dt  dt  dx ¥„ 


dt 

d 

dx//a 


d[(«;  1 

co 

1 

1  1 
o? 

F 

\L 

f 1  S  dcpa 

dx,  dy/a 

[\ 

v  p  dxt  dxt 

-  (Sa 

w)  pA  8 

(Scmp  1  F  P,  1 

''  ^  d¥a 

1  n 

A 

) 


x¥)  P, 


where, 


(58) 


SR  -d> 

a  a 


and 


Scmp  =  0 

a 


a  =  1 


SR  =-Y(h°W  co  )  and  Sc,np  =  ~(^  +  u ,^-  +  a  ^) 
°  ^  “  p  dt  'dx.  J  dx. 

>  i  i 


a  -  N 


5+1 


a= 1 


(59) 
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Equation  (58)  is  an  exact  transport  equation  for  the  scalar  FMDF  in  compressible  flows.  In  this 
equation,  the  Lewis  number  is  assumed  to  be  unity,  so  the  mass  and  thermal  diffusion 

coefficients  will  be  similarly  obtained  from  the  viscosity  as  T  =  /u  /  Sc  .  The  FMDF  equation 
cannot  be  solved  directly  due  to  the  presence  of  three  unclosed  terms.  Following  the  suggested 
models  for  these  unclosed  terms  [16-21],  the  closed  form  of  FMDF  transport  equation  for  a 
compressible  reacting  system  is  achieved, 


dPL  ,  d[{ui)LPL]=d_ 

dx, 


dt 


-  +  - 


dx, 


(r+r,) 


o(PLl{p) 


dx, 


(60) 


where. 


SR  =  co 

a  a 


and 


SR  =  -Y(h°W  a  )  and  Scmp  = 

a  v  a  a  a'  a 


a= 1 


a 

Hi  dt  dx. 


+  \T.  +(  <7 


,J  \  j/l]  dx 


d(U<)L) 


a  =  1 


a-N 


(61) 


where r,  =  (p'j /  vr  / Pr(  is  the  turbulent  diffusivity  and  Prt  is  the  turbulent  Prandtl  number 

(turbulent  Prandtl  and  Schmidt  numbers  are  the  same).  The  SGS  mixing  frequency  £lm  is 
evaluated  using  the  following  relation: 


(r+r,) 

<a!H) 


(62) 


Using  Equation  (60),  the  equation  for  the  first  SGS  moment  ((</>) L  )  is  derived  in  the  following 
form: 


_d_ 

dx. 

i 


f 

[(r  +  r) 


+  (p)  scmp 

Y  /  l  a 


(63) 


This  equation  can  also  be  obtained  directly  by  filtering  Equation  (57),  using  a  standard  gradient 
model  for  the  subgrid  flux  terms  and  neglecting  the  SGS  viscous  and  pressure  terms. 
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3.3.3  Hydrogen-Air  Reaction  and  Thermophysical  Properties 

Basic  Definitions 

Multi-step  reactions  are  usually  defined  by  the  following  equations: 


1  =  1,2,- N 


(64) 


where  M;  is  the  chemical  symbol  for  the  ith  species  and  v;,/  denotes  the  molar  concentration 
coefficient  of  the  ilh  species  in  reaction  I.  The  forward  and  backward  rate  constants  kf  and  kb  can 
be  calculated  using  the  Arrhenius  law,  as  follows: 


k(T)  =  AT"  exp (~Ea  /  (R°T)) 


(65) 


The  rate  of  a  reaction  depends  on  the  rate  constants  and  the  concentrations  of  species  involving 
in  the  reaction. 


co,  =  k, 


Ns 

no’ 


A  ||< 


1  =  1,2,-,  AT, 


i=l 


i= 1 


(66) 


The  net  rate  of  change  of  concentration  of  a  species  in  the  whole  reaction  set  is  calculated  from 
the  following  relation: 


dC  ^ 

&  =  —  =  I>"  ~  Ki>i 

i=i 


dt 


i  =  l,2,-,N 


(67) 


The  concentration  and  mass  fraction  of  the  i,h  species  are  related  through  the  following  equation: 

(68) 


C,  =  (i)  1  P 

1  v  TTZ  7 


W‘  f  (A  RT 


Subscript  “mix”  stands  for  the  mixture. 

Reaction  Mechanisms 

There  are  many  reaction  mechanisms  for  the  hydrogen-air  combustion.  Four  of  these 
mechanisms  are  considered  in  the  present  study.  These  mechanisms  are  presented  in  Tables  1 
through  4.  The  37-step  detailed  [8]  and  the  16-step  skeletal  [9]  mechanisms  solve  for  nine 
species,  while  the  14-step  skeletal  [10]  and  the  4-step  reduced  [9]  mechanisms  include 
contributions  of  eight  and  seven  species,  respectively.  The  species  H2O2  is  assumed  to  be  in 
steady  state  condition  in  the  last  two  mechanisms.  In  the  4-step  reduced  mechanism,  the  steady- 
state  assumption  is  made  for  HO2  as  well.  In  some  reactions  a  third  body,  M,  participates  in  the 
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reaction.  The  third  body  is  needed  to  carry  away  the  extra  energy  and  remains  unchanged  on 
both  reactant  and  product  sides  of  the  reaction.  In  all  reactions,  Mhas  a  fixed  efficiency  given 

by, 

M  =  1.0  H2  +  6.5  H20  +  0.4  02  +0.4  N2  +  1.0  O  +  1.0  H  +  1.0  OH  +  1.0  H02  +  1.0  H202  (69) 

It  must  be  noted  that  the  only  effect  of  variation  in  the  efficiencies  of  the  third  body  is  changing 
the  ignition  delay. 


Table  1:  The  37-Step  Detailed  Mechanism 


# 

Reaction 

# 

Reaction 

# 

Reaction 

1 

02  +  H  — »  OH  +  O 

14 

02  +  M— >20  +  M 

27 

H02  +  H02  — *  h2o2  +  o2 

2 

OH  +  O  — >  02  +  H 

15 

H  +  02  +  M  ->•  H02  +  M 

28 

2  OH  +  M  — >  H202  +  M 

3 

H2  +  O  — >  OH  +  H 

16 

ho2  +  m^h  +  o2  +  m 

29 

H202  +  M  — >  2  OH  +  M 

4 

OH  +  H  — >  H2  +  O 

17 

HO,  +  H  -►  OH  +  OH 

30 

H202  +  H  -*  H2  +  H02 

5 

H2  +  OH  — >  H20  +  H 

18 

OH  +  OH  — >  H02  +  H 

31 

H2  +  H02  ->  H202  +  H 

6 

H20  +  H  ->  H2  +  OH 

19 

H02  +  H  ->  H2  +  02 

32 

H202  +  H  -►  H20  +  OH 

7 

OH  +  OH  — »■  H20  +  O 

20 

H2  +  02  H02  +  H 

33 

H20  +  OH  ->  H202  +  H 

8 

H20  +  O  ->•  OH  +  OH 

21 

ho2+h->h2o  +  o 

34 

H202  +  O  OH  +  ho2 

9 

2H  +  M— >H2  +  M 

22 

H20  +  O  ->  H02  +  H 

35 

OH  +  H02  H202  +  O 

10 

H2  +  M— »2H  +  M 

23 

H02  +  O  ->•  OH  +  02 

36 

H202  +  OH  ->  H20  +  H02 

11 

H  +  OH  +  M  ->  H20  +  M 

24 

OH  +  02  ->  H02  +  O 

37 

H20  +  H02  ->  H202  +  OH 

12 

H20  +  M^H  +  OH  +  M 

25 

H02  +  OH  ->  H20  +  02 

13 

20  +  M— >02  +  M 

26 

H20  +  02  -►  H02  +  OH 

21 
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Table  2:  The  16-Step  Skeletal  Mechanism 


# 

Reaction 

# 

Reaction 

# 

Reaction 

1 

H  +  02  — >  O  +  OH 

7 

h  +  h  +  m->h2  +  m 

13 

H02  +  H  ->  OH  +  OH 

2 

O  +  OH  — >  H  +  02 

8 

H  +  O  +  M  — >  OH  +  M 

14 

ho2  +  ho2  — *  h2o2  +  o2 

3 

O  +  H2  — >  H  +  OH 

9 

H  +  OH  +  M  — >  H20  +  M 

15 

H202  +  M  ->  OH  +  OH  +  M 

4 

H  +  OH  — >  O  +  H2 

10 

H  +  02  +  M  — >  H02  +  M 

16 

H02  +  H2  H202  +  H 

5 

OH  +  H2  — >  H  +  H20 

11 

H02  +  H  -*  H2  +  02 

6 

H  +  H20  -►  OH  +  H2 

12 

H2  +  02  — ^  H02  +  H 

Table  3:  The  14-Step  Skeletal  Mechanism 


# 

Reaction 

# 

Reaction 

# 

Reaction 

1 

H02  +  H  -»  H2  +  02 

6 

H2  +  OH  — ►  H20  +  H 

11 

H  +  OH  +  M  -*  H20  +  M 

2 

OH  +  O  — >  02  +  H 

7 

H20  +  H  ->  H2  +  OH 

12 

H  +  02  +  M  — >  H02  +  M 

3 

02  +  H  — >  OH  +  O 

8 

OH  +  OH  — >  H20  +  O 

13 

H02  +  OH  ->•  H20  +  02 

4 

H2  +  O  — »  OH  +  H 

9 

H20  +  O  -»  OH  +  OH 

14 

H2  +  02  H02  +  H 

5 

OH  +  H  — >  H2  +  O 

10 

H02  +  H  OH  +  OH 

Table  4:  The  4-Step  Reduced  Mechanism 


# 

Reaction 

# 

Reaction 

1 

H2  H  +  H 

3 

H2  +  OH  <-»  H20  +  H 

2 

H  +  02  0  +  OH 

4 

OH  +  OH  <->  H2  +  02 

Thermodynamic  and  Transport  Properties 

In  this  work,  the  constant-pressure  heat  capacity  (Cp)  of  the  ith  species  is  evaluated  using  the 
NASA  polynomial  in  the  following  form: 


~r  ~  a\j  +  a2 +  a3,;^2  + 

K 


(70) 
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This  method  is  also  used  in  GRI  Mech  and  CHEMKIN-II  codes.  The  produced  data  from  the 
above  formula  are  read  from  an  input  in  the  main  code.  By  keeping  the  temperature  intervals 
constant,  there  will  be  a  one-to-one  correspondence  of  temperature  and  the  value  of  Cp  for  each 
component.  This  method  is  very  efficient  and  reduces  the  run  time. 

When  the  Cp  of  each  species  is  calculated,  the  mixture  Cp  can  be  obtained  as  follows: 

Ns 

(71) 

7=1 

There  are  three  different  methods  for  calculating  the  molecular  viscosity  and  thermal 
conductivity  of  the  species.  In  the  first  method,  one  can  calculate  the  viscosity  and  conductivity 
of  the  ith  species  using  the  following  formulas  proposed  by  NASA  [11]: 

ln(M.)  =  41nr  +  ^  +  ^  +  JDft  (72) 

Ba  Cx  (73) 

ln(Ai)  =  A,\nT  +  ^  +  ^  +  DAi 


The  second  method  proposes  two  theoretical  relations  based  on  the  gas  kinetic  theory  to  calculate 
the  viscosity  and  conductivity  of  the  species  [9]. 


A, 


8.44107xl(T8 


JWT 


(74) 


L  = 


--+C, 

4  m; 


A, 


(75) 


In  these  relations,  Q12'2'  is  the  collision  integral  and  depends  on  the  nondimensional  temperature 
T  and  parameter  5  .  In  addition,  k  ,  a,,  and  in,  are  Boltzmann  constant  (=1.3806x10  J/K), 

Lennard-Jonnes  collision  diameter  and  mass  of  the  ith  species  within  the  mixture. 

In  the  third  method,  one  can  fit  polynomials  to  the  data  of  CHEMKIN  to  obtain  the  following 
relations  to  calculate  viscosity  and  thermal  conductivity  of  the  ith  individual  component: 


In(A;)  =  IX7  [ln(r)] 

7=0 

(76) 

Npo,  J 

ln«)=2>,v[lnm] 

(77) 

7=0 


Where,  Npoi  is  the  order  of  the  polynomial. 
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The  other  important  transport  property  is  the  binary  diffusion  coefficient  between  species  i  and  j 
which  is  given  by  the  following  theoretical  relation  in  the  standard  units: 


D..  =5. 87646x10  6 


Jt3(wi.+w/)/wiw/ 


(78) 


where  Q(1,1>  is  the  collision  integral,/;  is  the  pressure  (in  tenns  of  atmosphere),  and  a,j,  is  the 
effective  Lennard-Jones  collision  diameter.  Since  Djj  is  only  a  function  of  temperature  and 
pressure,  we  can  fit  polynomials  to  the  (Djj  .  p)  and  then  use  them  in  the  main  code.  This  will 
result  in  a  significant  reduction  of  computational  cost  without  loss  of  accuracy.  A  third-order 
polynomial  in  the  following  format  is  found  to  be  quite  satisfactory: 

ln(Z>.  •  p)  =  d0  ij  +  dUJ  HT)  +  d2  ij  [ln(T)f  +  d3  ij  [ln(T)]'  (79) 


There  are  different  ways  of  calculating  the  transport  properties  of  the  mixture  from  the  properties 
of  the  individual  components.  The  exact  theoretical  approach  yields  the  following  relations  to 
calculate  mixture  viscosity  and  conductivity  [9]: 


(80) 

(81) 


(82) 


Since  there  are  two  exterior  and  interior  summations  in  these  relations,  their  computational  costs 
are  relatively  high.  Thus,  we  look  for  a  suitable  approximate  approach  that  leads  to  reasonable 
results  and  is  computationally  cheaper  at  the  same  time.  Three  different  approximations  are 
separately  listed  below  for  the  mixture  viscosity  and  conductivity. 

First  approximation  for  pmix  (proposed  by  Wilke  [12]  and  modified  by  Bird,  et  al.  [13]): 


(=i 


XjPi 


j=i 


(83) 
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Second  approximation  for  jumix\ 


Ns 

li  .  nfjT.it.  (84) 

•  mix  ir  i  v  ' 

i=l 

Third  approximation  for  jumix: 

Ns 

fu.UyYfu  (85) 

r  mix  ir  i  x  ' 

i=l 

Similarly,  we  have  the  following  approximations  for  the  mixture  conductivity: 

First  approximation  for  Xmix  (proposed  by  Mathur,  et  al.  [14]): 


= 


2 


(86) 


Second  Approximation  for  Xmix: 

Ns 

2.0V/.1.  (87) 

mix  i  i  v  ' 

i= 1 

Third  Approximation  for  Xmix: 

Ns 

JY  Ya.  (88) 

mix  i  i  v  y 

/= 1 

3.4  Implementation  of  FMDF  in  the  DSEM  Code 

This  section  is  dedicated  to  the  implementation  of  the  FMDF  method  in  the  two-dimensional 
DSEM  code.  Previously,  the  DSEM  code  contained  routines  for  tracking  inertial  (heavy) 
particles;  however,  to  be  compatible  with  the  FMDF  method,  it  is  necessary  to  modify  these 
existing  routines  and  add  new  routines  to  account  for  additional  terms  and  equations.  The  plan 
for  implementation  is  as  follows: 

•  Modeling  of  Monte  Carlo  particles 

•  Validation  of  particle  routines 

•  Development  of  FMDF  routines  compatible  with  DSEM  code 

•  Calculation  of  Eulerian  field  from  tracked  particles 

•  Consistency  test  comparing  DSEM  results  to  FMDF  results 
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3.4.1  Mathematical  Modeling 

In  order  to  explain  the  implementation  of  the  particle  tracking  and  FMDF  algorithms,  we  first 
present  the  necessary  stochastic  differential  equations.  For  a  detailed  explanation  of  how  particle 
tracking  is  implemented  in  the  DSEM  code  we  refer  to  [23].  For  details  of  the  derivation  of  the 
FMDF  equations  we  refer  to  [17]. 

Equation  of  Motion  fora  Particle  without  Inertia 

The  Lagrangian  procedure  used  to  implement  the  tracking  of  particles  in  the  DSEM  code  is  a 
modification  of  the  method  described  by  [23],  Considering  that  the  FMDF  method  is  concerned 
with  particles  with  no  inertia,  the  governing  equations  of  motion  for  our  particles  are  simplified 
to  a  single  equation  for  particle  position 


dXp 

dt 


=  vf  +  Dt 


(89) 


where  Xp  is  the  position  of  the  particle  in  space,  vp  is  the  velocity  of  the  carrier  phase  at  the 
particle  location,  and  Dt  accounts  for  the  diffusion  process.  Such  particles  do  not  have  any  mass 
and  therefore  should  behave  as  fluid  particles  of  the  carrier  phase  rather  than  inertial  particles 
entrained  in  the  flow.  It  should  be  noted,  however,  that  particles  with  mass  and  inertia  behave 
differently,  for  this  purpose  a  differential  comparison  is  performed  to  ensure  that  the  particles  we 
are  modeling  behave  as  expected.  The  complete  set  of  Lagrangian  equations  for  inertial  particles 
are  found  in  [23]. 

Additional  Stochastic  Differential  Equations 

In  FMDF,  the  compositional  values  of  the  particles  are  changed  due  to  their  mixing  and  reaction 
within  the  flow.  For  any  given  scalar  <p  the  subgrid  scale  reaction  may  be  described  as 


d<fia 

dt 


(90) 


whereas  the  subgrid  scale  mixing  is 


d(f>a 

dt 


^ m(0a  (0a)) 


(91) 


where  0a=0a(Xp(t),  t)  is  a  scalar  attributed  to  the  particle  at  the  particle  location,  Xp(t),  in 
space,  and  <  >  shows  the  ensemble  average.  In  addition,  the  term  S(p  describes  the  reaction 
source,  and  nmis  a  modeled  tenn  relating  to  the  frequency  of  mixing  at  the  subgrid  level. 

3.4.2  Numerical  Implementation 

To  successfully  implement  the  FMDF  method  into  the  DSEM  code,  it  is  first  necessary  to 
develop  a  particle-tracking  algorithm  suitable  to  the  DSEM  method.  The  algorithm  must  contain 
the  following  steps: 

•  Seed  particles  by  initializing  their  position  and  velocity 

•  Use  search  algorithm  to  locate  particles  within  elements 
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•  Map  particles  from  Chebyshev  grid  to  an  equidistant  grid 

•  Find  carrier  phase  properties  at  particle  location  using  interpolation  scheme 

•  Integrate  particle  equations  in  time 

•  Establish  ensemble  domains  around  Eulerian  grid  points 

•  Take  statistics  from  ensemble  domains 

•  Calculate  Eulerian  fields  based  on  Lagrangian  fields 

In  this  section,  we  briefly  present  the  necessary  work  performed  to  implement  the 
aforementioned  steps  into  the  DSEM  code. 

Search  Algorithm 

Considering  the  algorithm  in  use  tracks  each  particle  individually  in  a  Lagrangian  manner,  it  is 
necessary  to  locate  individual  particles  within  our  elements.  Locating  the  particles  is 
accomplished  by  scanning  the  entire  computational  domain.  For  the  case  of  subdomains  with 
rectangular  shapes,  such  as  in  the  shear  flow  case  presented  below;  it  is  possible  to  take 
advantage  of  an  orthogonal  grid  in  mapped  space  to  reduce  computational  cost.  Once  the  particle 
location  is  detennined  in  mapped  space  on  the  Chebyshev  grid,  it  must  then  be  mapped  to  an 
equidistant  grid  where  the  particle  equations  shall  be  integrated  in  physical  space. 

Interpolation  Scheme 

The  DSEM  method  provides  solution  values  at  Gauss-Gauss  points.  However,  it  becomes 
instantly  evident  that  the  particles,  which  we  are  tracking,  do  not  lie  upon  these  points.  It  is 
therefore  necessary  to  evaluate  the  properties  of  the  carrier  phase  at  the  particle  location  in  space. 
Although  there  are  several  methods  available  to  perform  such  operations,  the  Lagrangian 
interpolation  scheme  of  order  six  is  selected  due  to  its  low  cost  and  high  accuracy.  [23] 

The  sixth  order  Lagrangian  interpolation  scheme  uses  values  taken  from  three  grid  points  on 
either  side  of  the  particle  location  in  mapped  space  to  calculate  a  polynomial.  A  value  at  the 
particle  location  is  then  calculated  from  the  resulting  sixth-order  polynomial.  If  the  particle  is 
located  near  the  boundary  of  the  element,  we  still  use  all  the  interpolation  points  from  the  same 
element.  We  have  shown  that  this  does  not  result  in  a  loss  of  accuracy  due  to  the  fact  that  the 
Chebyshev  points  are  distributed  more  closely  to  each  other  near  the  element  boundary.  The  fact 
that  the  interpolation  scheme  is  based  entirely  on  the  points  within  the  same  element  significantly 
reduces  the  cost  of  message  passing  in  MPI  and  thus  adds  to  the  efficiency  of  parallelization. 

Time  Stepping  Scheme 

The  temporal  discretization  of  the  transport  equation  is  done  using  an  Adams-Bashforth  (AB) 
scheme.  In  order  to  employ  the  AB  scheme,  all  variables  for  the  particle  properties  must  be 
projected  onto  the  physical  space  from  the  mapped  space.  The  interpolation  scheme  provides  us 
with  the  necessary  information  to  solve  the  velocity  at  the  particles  location  in  physical  space 

(92) 


(93) 
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where  Vf  —  (vx,  vy  )is  the  fluid  velocity  at  the  particle  location  and  Q1  and  Q2  are  defined  in 
Equation  (30).  To  find  the  velocity  of  the  massless  particle  in  the  mapped  space,  we  simply  use 

Vp  =  vf  (94) 

where  n  signifies  the  time  level.  Knowing  the  velocity,  we  then  update  the  particle  position  as 


dXp 

dt 


=  Vf,  X 


n+ 1 
P 


—  Xp  +  ^  At(3i7p 


(95) 


Ensemble  Domain 

To  solve  the  FMDF  equations,  a  sampling  domain  is  constructed  from  which  we  are  able  to  take 
statistics.  Such  a  domain  is  referred  to  as  an  “ensemble  domain.”  Ideally  the  ensemble  domain 
would  be  infinitely  small,  and  the  number  of  particles  within  the  ensemble  domain  would  be 
infinitely  large.  However,  computationally  we  are  restricted  to  use  a  finite  size  in  space  and  a 
finite  number  of  particles.  The  challenge  is  to  find  the  smallest  number  of  particles  that  would 
provide  ensemble-averaged  quantities  without  a  significant  error. 

The  solutions  for  DSEM  and  FMDF  are  treated  in  a  split  manner.  The  FMDF  method  only 
carries  dependency  on  the  velocity  field  calculated  in  the  DSEM  and  the  primitive  variables 
assigned  to  the  particles  when  they  are  initialized.  However,  due  to  the  fundamental  differences 
between  the  DSEM  method  and  LES-FD  method,  as  shown  in  [17],  it  was  necessary  to  create  an 
additional  Lagrangian  sub-grid  on  which  the  FMDF  solutions  would  lie.  For  the  purposes  of  this 
consistency  test  a  uniform  grid  was  selected  as  shown  in  Figure  4,  but  a  nonuniform  grid  may  be 
easily  substituted  if  necessary. 
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Figure  3:  DSEM  Grid  for  a  Free  Shear  Flow 


Figure  4:  Ensemble  Grid  for  a  Free  Shear  Flow 
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Initial  and  Boundary  Conditions 

The  particles  are  initialized,  or  rather  seeded,  consistently  with  respect  to  the  initial  flow  field, 
meaning  they  carry  the  initial  velocity,  temperature,  and  weighting  as  the  fluid  at  their  location  in 
mapped  space.  The  particles  should  be  distributed  uniformly,  yet  randomly,  within  the  domain  of 
the  carrier  phase. 

In  addition,  should  particles  leave  the  flow  field  during  simulation,  they  are  replaced  at  the  inlet 
to  ensure  that  a  consistent  number  of  particles  remain  within  the  domain  for  averaging. 

Weighting  becomes  necessary  with  density  variation  in  the  carrier  phase  as  particles  may  be 
sparser  in  such  areas.  To  ensure  proper  averaging  occurs,  additional  particles  shall  be  added  to 
areas  of  lower  density,  and  shall  be  weighted  accordingly. 

In  order  to  reduce  computational  cost  and  memory  consumption,  a  method  of  variable  weighting 
is  employed.  Considering  that  in  the  FMDF  method  each  particle  is  representative  of  a  fluid 
element  within  the  domain,  we  are  able  to  use  a  system  of  weights  that  are  proportional  to  the 
mass  of  the  element  in  question.  Doing  so  allows  for  a  smaller  number  of  particles  to  be  tracked. 
Computation  of  the  weights  of  the  aforementioned  particles  is  trivial  and  is  dependent  on  the 
initial  number  of  particles  populating  an  ensemble  element  as  well  as  the  fluid  density  at  that 
specific  location  in  space.  An  individual  particle’s  mass  is  calculated  as  an  initial  condition  as 
such: 


m 


mp  = 


n 


ae 


(96) 


where  npE  is  the  number  of  particles  contained  within  the  ensemble  domain  and  m^E  is  the  mass 
of  fluid  contained  within  the  ensemble  element  encasing  the  particle  in  question.  The  mass  of  the 
fluid  is  simply  found  by  the  relation  m^E  =  PfV-,  where  Pj  is  the  fluid  density  and  V  is  the 
element  volume.  It  should  be  emphasized  that  the  particle  mass  is  used  for  keeping  track  of 
density  only.  For  the  purpose  of  tracking  the  trajectories,  all  particles  are  considered  massless. 

To  assign  initial  properties  we  take  advantage  of  the  DSEM  code’s  Lagrangian  interpolation 
routines  to  assign  velocity,  density,  and  a  scalar  to  the  particle  from  the  Eulerian  field.  However, 
after  the  initial  conditions  are  assigned,  only  the  velocity  continues  to  be  interpolated  onto  the 
particle  location.  The  weighting  is  to  remain  the  same  unless  additional  particles  need  to  be 
added  or  particles  need  to  be  removed,  and  the  scalar  is  to  be  solved  by  the  FMDF  stochastic 
differential  equations.  It  can  clearly  be  seen  that  the  two  methods,  although  being  used  together, 
are  in  fact  providing  independent  solutions. 

To  achieve  appropriate  comparison,  we  require  that  the  boundary  conditions  imposed  on  both 
solvers  be  equal.  This  is  done  by  matching  the  mass  flow  rate  of  fluid  in  the  Eulerian  field  with 
the  mass  flow  rate  of  particles  entering  at  the  inlet: 


rhiniet  — 


(97) 


which  in  the  appropriate  variables,  when  solved  for  number  of  particles  to  be  seeded  at  a  given 
time  step,  reduces  to 
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PfUfAAeAt 


(98) 


Fluid  density,  Pf,  and  fluid  velocity,  Uf,  are  taken  at  the  location  of  the  ensemble  into  which  the 
particles  are  to  be  seeded.  The  inlet  area,  AAr,  is  taken  as  the  area  (or  for  the  case  of  two- 
dimensions  it  is  the  length)  of  the  ensemble  element.  Since  the  number  of  particles  being  added 
may  only  be  an  integer  number,  there  are  two  available  options  for  handling  the  remainder.  One 
may  either  adjust  the  weighting  of  the  particles  to  match  the  mass  per  time  step  or  perform 
integer  subtraction  and  track  remainders  which  are  carried  over  to  future  time  steps.  To  minimize 
computational  overhead  from  additional  particles  with  lower  weighting,  integer  additions  with 
fixed  weighting  were  selected.  Although  this  causes  the  mass  flow  rates  at  each  individual  time 
step  to  be  unequal,  the  time  average  of  mass  entering  is  equal  and  the  method  holds  true. 

Calculation  of  Scalar  Density 

Due  to  the  fact  that  the  particles  represent  discrete  fluid  elements  in  our  simulations,  it  is 
expected  that  if  fluid  elements  enter  or  leave  a  respective  ensemble  domain  the  density  of 
aforementioned  domain  will  increase  or  decrease.  This  is  shown  by  the  FMDF  formulation  [17] 
as: 


<P>  =  C  ^  Wn 

n  E  Ae 


(99) 


where  C  is  a  proportionality  constant  and  Wn  is  the  weighting  of  particle  n  that  is  within  an 
element  of  the  ensemble  domain,  AE.  It  is  seen  that  the  density  of  an  ensemble  domain  is  in  fact 
proportional  to  the  sum  of  the  weights  of  the  Monte  Carlo  particles  being  tracked.  The  constant 
may  be  found  a  posteriori.  For  the  case  of  unifonn  weights,  this  equation  reduces  to  (p)  ~ 

—  np.By  using  uniform  weights,  as  density  decreases  in  an  ensemble  volume  the  particle 

quantity  shall  also  decrease  which  in  turn  decreases  statistical  accuracy.  It  is  for  this  purpose  that 
particle  adding/removal  with  variable  weighting  is  implemented  in  the  DSEM/FMDF  code. 
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4.0  Results  and  Discussions 

4.1  Entropy  Viscosity  Method  in  1 D 

As  mentioned  earlier,  we  consider  two  model  problems  in  this  work  to  validate  our  code  and  see 
how  the  entropy  viscosity  method  and  explicit  filtering  can  be  effective  in  removing  spurious 
oscillations  and  shocks.  Linear  advection  and  nonlinear  inviscid  Burgers  equations  are 
considered  here  and  they  are  solved  with  different  initial  conditions  ranging  from  smooth  waves 
to  extremely  sharp  waves  to  observe  how  the  code  handles  different  types  of  initial  conditions. 

4.1.1  One-Dimensional  Linear  Advection  Equation 

The  model  problem  is 

ut  +  ux  =  0  in  fl  =  [0,1]  (100) 

with  periodic  boundary  condition.  Three  different  types  of  initial  conditions  are  considered  in 
this  study,  namely  sine  wave,  Gaussian  exponential  wave  and  square  wave.  We  consider  these 
three  waves  as  three  different  levels  of  smoothness  of  initial  disturbance.  The  simulations  have 
been  conducted  within  a  domain  with  twenty  P8  elements,  where  P  indicates  the  order  of  the 
approximation  polynomial  within  each  element.  The  results  are  depicted  in  Figure  5  for  different 
types  of  initial  conditions.  In  Figure  5(a),  the  solution  without  using  filtering  or  entropy  viscosity 
is  depicted  in  comparison  with  the  exact  solution.  As  it  is  clearly  shown,  for  the  sine  wave  the 
unaltered  DSEM  method  is  able  to  generate  a  solution  with  a  good  agreement  to  exact  solution. 
For  this  smooth  initial  condition,  there  are  no  oscillations  in  the  results  after  20  time  units.  We 
can  also  use  this  case  to  assess  the  numerical  diffusion  introduced  by  the  filtering  and  the  EV 
method  using  the  results  displayed  in  Figure  5(b), (c).  By  filtering  from  P8  to  P5,  we  cannot  see 
any  significant  diffusion  error,  although  a  slight  dispersion  error  is  visible  after  20  periods,  as 
evident  by  the  slight  shifting  of  the  wave  to  the  right  of  the  exact  solution.  In  contrast,  applying 
EV  damps  the  amplitude  of  the  wave  by  about  5%  and  this  may  be  interpreted  as  a  measure  of 
the  diffusion  error.  There  is  also  a  visible  dispersion  error  associated  with  EV,  larger  than  that 
introduced  by  the  filtering. 

As  the  next  test  case,  we  have  combined  a  Gaussian  exponential  wave  and  a  square  wave  to 
analyze  the  stability  and  accuracy  of  the  DSEM  method  in  comparison  to  the  results  obtained 
after  using  filtering  or  EV.  Figure  6(a)  shows  a  comparison  between  the  exact  solution  and  the 
unaltered  DSEM  results.  As  seen  from  the  figure,  after  5  time  units  relatively  large  oscillations 
appear  on  top  and  around  the  square.  It  is  clear  that  the  unaltered  DSEM  method  is  not  capable 
of  handling  abrupt  discontinuities  in  the  solution  and  we  need  some  remedy  to  overcome  this 
problem.  In  Figure  6(b)  and  2(c)  the  results  after  applying  filter  and  EV  are  shown.  While 
filtering  removes  the  high-frequency  oscillations,  it  is  clearly  ineffective  in  providing  a 
completely  oscillation-free  solution.  However,  we  can  clearly  see  that  EV  is  able  to  damp  the 
oscillations  and  smoothen  the  solution  without  incurring  a  significant  loss  of  accuracy. 
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(c)  DSEM  with  Entropy  Viscosity 


Figure  5:  Solution  of  the  Linear  Advection  Equation  for  an  Initial  Sine  Wave  in  a  Domain  with  20  P8 

Elements  at  t=20 


Figure  6:  Solution  of  the  Linear  Advection  Equation  for  an  Initial  Combined  Square  and  Exponential  Wave 

in  a  Domain  with  20  P8  Elements  at  t=5 
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4.1.2  One-Dimensional  Inviscid  Burgers  Equation 

The  model  problem  for  this  case  is 

ut  +  uux  —  0  in  H  =  [0,1]  (101) 

with  periodic  boundary  conditions.  We  study  this  case  for  three  different  initial  conditions  and 
compare  the  solutions  with  and  without  EV.  Filtering  proved  not  to  be  helpful  in  removing 
oscillations  in  this  case,  as  we  also  saw  in  linear  advection  problem.  Figure  7(a)  shows  the  results 
for  a  sine  wave  initial  condition  i/0(x)  =  sin(27rx).  As  seen  in  the  figure,  the  initial  sine  wave  is 
being  deformed  in  time  until  a  shock  forms  at  x  —  0.5.  The  unaltered  DSEM  code  is  capable  of 
handling  this  problem  smoothly  with  no  oscillations.  Results  from  EV  method  are  also  stable 
with  negligible  deviations  from  the  exact  solution  because  of  the  added  diffusion. 


t 

0.5 

0 

5 


-0.5 


-1 


(a)  Unaltered  DSEM 


(b)  DSEM  with  Entropy  Viscosity 


Figure  7:  Solution  of  the  Inviscid  Burgers  Equation  for  an  Initial  Sine  Wave  with  20  P16  Elements 


Figure  8:  Solution  of  the  Inviscid  Burgers  Equation  for  an  Initial  Exponential  Wave  with  20  P16  Elements 


In  Figure  8  results  are  presented  for  the  following  initial  data  with  and  without  using  the  EV 
method: 


Uq 


_|_  g  — 300(x— 0.3)2 


if  |x  —  0.3 1  <  0.15 
otherwise 


(102) 
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The  figure  shows  that  the  wave  initially  starts  tilting  to  the  right  and  as  soon  as  the  shock  forms  it 
starts  moving.  It  is  obvious  that  the  DSEM  method  is  not  able  to  capture  the  dynamics  of  this 
shock  and  produces  severe  oscillations.  In  contrast,  Figure  8(b)  shows  that  enhancing  DSEM 
with  EV  substantially  reduces  the  magnitude  of  the  oscillations  and,  with  the  exception  of  some 
small  wiggles  near  the  shock,  the  solution  can  otherwise  be  considered  smooth.  It  is  also 
worthwhile  to  mention  that  the  remaining  oscillations  would  not  grow  in  time  in  the  case  with 
EV,  while  the  big  oscillations  in  Figure  8(a)  grow  in  time  and  we  are  not  able  to  obtain  a  solution 
for  longer  times. 
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Figure  9:  Solution  of  the  Inviscid  Burgers  Equation  for  an  Initial  Exponential  Wave  with  20  P16  Elements 


As  the  last  case  we  consider  a  square  wave  as  initial  condition  and  this  is  the  most  abrupt  profile 
we  are  going  to  test  the  method  with.  As  we  can  see  in  Figure  9(a),  even  for  a  short  time  (t=0. 15) 
the  unaltered  DSEM  method  produces  very  oscillatory  solution  and  is  not  stable.  Figure  9(b) 
shows  the  results  for  DSEM  enhanced  with  EV  for  three  different  times  and  it  is  clearly  seen  that 
the  solution  is  generally  smooth  and  stable;  the  small  wiggles  behind  the  shock  are  not  growing 
in  time.  We  believe  that  these  remaining  oscillations  can  be  removed  with  further  fine-tuning  of 
the  method.  It  is  also  important  to  mention  that  the  expansion  wave,  on  the  left,  is  simulated 
accurately  with  or  without  EV. 

4.2  Entropy  Viscosity  Method  in  2D 

To  test  the  perfonnance  and  accuracy  of  the  developed  2D,  DSEM-EV  code  a  moving  shock 
problem  is  considered.  Computations  show  that  the  DSEM-EV  method  is  able  to  capture  shocks 
while  removing  spurious  Gibbs  oscillations.  Through  a  parametric  study  into  the  effect  of  grid 
resolution,  we  show  that  shocks  are  captured  sharply  with  the  EV  method  within  only  a  few 
quadrature  points,  requiring  less  resolution  than  lower-order  Roe  solver. 

Treatment  of  the  boundaries  and  interfaces  plays  a  major  role  in  obtaining  a  stable  and  accurate 
solution.  Figure  10  shows  a  comparison  between  the  solution  obtained  using  the  unaltered 
DSEM  Euler  solver  and  the  solution  obtained  from  EV-DSEM  Euler  solver  without  having 
boundary  and  interface  treatment.  Although  the  oscillations  tend  to  disappear  after  using  EV 
method  but  we  could  not  obtain  a  stable  solution  after  the  shock  reaches  the  interface.  After 
updating  the  EV-DSEM  code  to  treat  the  boundary  and  interfaces  this  issue  was  resolved. 
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Figure  10:  Velocity  and  Temperature  Profiles  with  20  Elements 


4.2.1  Moving  Shock  Problem 

As  we  mentioned,  updating  the  code  and  adding  the  feature  to  treat  the  boundary  and  interfaces 
to  enforce  the  C°  continuity  of  entropy  fluxes  we  were  able  to  achieve  stable  solution  for  the 
moving  shock  problem  using  EV-DSEM  Euler  solver.  To  demonstrate  this,  we  consider  a 
moving  shock  problem  in  2D.  Figure  1 1  shows  the  initial  and  boundary  conditions  as  well  as  a 
sample  grid.  A  slip  boundary  condition  is  applied  to  lower  and  upper  walls. 
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Figure  11:  Grid  and  Initial  Mach  Number  Distribution  for  the  Moving  Shock  Problem 


Figure  12:  Velocity  and  Temperature  Plots  for  the  Moving  Shock  Problem  with  2  P17  Elements  at  Different 

Times 
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As  seen  in  Figure  12,  a  smooth  decrease  in  velocity  field  is  considered  initially  at  x=l.l.  This 
velocity  field  transitions  into  a  shock  as  it  moves  towards  the  end  of  the  domain.  Figure  12 
shows  the  results  of  the  computation  with  only  two  elements  in  x-direction  and  a  very  high-order 
polynomial  of  17. 

It  is  clearly  seen  that  the  Gibbs  oscillations  are  entirely  suppressed,  and  the  shock  smoothly 
passes  through  the  interface  at  x=2.5.  Most  importantly,  the  solution  is  stable  and  smooth, 
whereas  computation  without  EV  is  unstable.  This  shows  a  proof-of-concept  of  the  EV  method. 
The  EV  method  dynamically  determines  an  appropriate  artificial  viscosity  at  the  shock  location, 
as  seen  in  Figure  13, in  regions  where  Gibbs  oscillations  are  present. 


Figure  13:  Sample  Plot  of  EV  Viscosity  and  its  Corresponding  Velocity  Field  for  a  Case  with  2  P17  Elements 

at  t  =  0.5 

It  is  of  interest  to  investigate  the  effect  of  influential  parameters  on  the  magnitude  of  the  error. 
We  first  investigate  the  effect  of  the  number  of  elements  in  the  domain. 

Effect  of  the  Number  of  Elements 

We  present  the  results  for  three  different  cases  at  a  certain  time.  We  consider  P4  elements  and 
show  the  results  for  grids  with  2,  4,  6  and  8  elements  in  the  domain.  Table  5  compares  the  RMS 
of  the  error  for  each  of  these  cases  relative  to  the  high-resolution  case  with  400  PI  elements, 
which  we  consider  as  the  reference  case  with  highest  accuracy.  As  we  can  see,  with  increasing 
the  number  of  elements  in  the  domain  the  error  decreases  drastically.  A  comparison  between  the 
velocity  profiles  obtained  with  the  EV  method  with  8  P4  elements  and  the  Euler  solver  (without 
EV)  with  400  PI  elements  is  shown  in  Figure  14.  This  comparison  reveals  a  good  agreement 
between  the  results  obtained  from  these  two  approaches.  As  expected,  the  application  of  the  EV 
method  reduces  the  sharpness  of  the  shock  slightly,  but  produces  a  stable  solution  with  a 
significantly  smaller  number  of  grid  points. 

Effect  of  the  Polynomial  Order 

A  similar  approach  is  used  to  study  the  effect  of  the  polynomial  order  on  the  error.  We  used  4 
elements  with  different  polynomial  orders  to  compare  the  error.  A  comparison  of  the  errors  in 
Table  5  and  Table  6  reveals  that,  for  the  same  total  number  of  grid  points,  using  higher  order 
polynomial  is  more  effective  in  decreasing  the  error  than  using  higher  number  of  elements.  This 
trend  is  expected  in  spectral  h/p  methods,  and  the  results  here  demonstrate  that  this  behavior  does 
not  change  with  implementation  of  the  EV  method. 
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3 


Figure  14:  Comparison  between  Velocity  Profiles  at  t  =  0.75  Obtained  Using  the  DSEM-EV  with  8  P4 
Elements  and  the  DSEM  (without  EV)  with  400  PI  Elements 

Table  5:  RMS  Error  with  P4  Approximation  for  a  Different  Number  of  Elements 


Number  of  Elements 

2 

4 

8 

16 

RMSE 

2.28E-1 

1.43E-1 

5.80E-2 

2.10E-2 

Table  6:  RMS  Error  Using  Four  Elements  with  Different  Polynomial  Orders 


Polynomial  order 

2 

3 

4 

6 

12 

RMSE 

2.82E-1 

1.81E-1 

1.04E-1 

5.10E-2 

1.61E-2 

4.3  Performance  of  Chemistry  Model 

To  compare  the  performance  of  the  different  mechanisms  and  models,  a  simple  reaction  problem 
is  considered.  This  problem  is  a  perfectly-stirred,  premixed,  iso-pressure  reactor.  The  species  and 
energy  equations  for  this  problem  are  given  by, 


%pYt) 

dt 


□ 

=  W.  mi 

l 


i  =  1,2, 


N 


(103) 


d(phs) 

dt 


NS  3 


V  h°  m 


(104) 


In  the  latter  equation,  hs  and  h  are  sensible  and  formation  enthalpies,  respectively.  The 
numerical  solution  begins  from  an  initial  temperature  and  species  mass  fractions  and  then  the 
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two  mentioned  equations  are  solved  explicitly  for  new  values  of  temperature  and  species  mass 
fractions.  As  long  as  the  time  step  is  kept  small  enough,  the  explicit  solver  can  be  used. 

The  results  obtained  for  different  mechanisms  are  compared  in  Figures  13-17.  According  to  these 
figures,  the  4-step  reduced  mechanism  acts  completely  different  from  the  other  mechanisms.  The 
two  skeletal  mechanisms  have  mostly  the  same  performances.  All  mechanisms  except  the  4-step 
one,  predict  the  same  value  for  the  ignition  delay. 


\ 

\ 
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Figure  15:  Temperature  vs.  Time  for  Different  Mechanisms,  (cp  =1.0) 


Figure  16:  Mixture  Density  vs.  Time  for  Different  Mechanisms,  (cp  =  1.0) 
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Figure  17:  Mole  Fraction  of  H2  vs.  Time,  (cp  =  1.0) 


Figure  19:  Mole  Fraction  of  OH  vs.  Time,  (<p  =  1.0) 


The  results  obtained  by  our  in-house  code  are  compared  to  those  of  CHEMKIN  in  this  section. 
Although  the  reaction  parameters  are  not  identical  for  the  compared  models,  we  have  found  the 
results  obtained  by  our  in-house  chemistry  code  to  be  the  same  as  those  of  the  CHEMKIN.  This 
is  illustrated  in  Figures  18-22,  where  the  temporal  variation  of  the  temperature,  mole  fractions  of 
some  species,  and  mixture  density  for  the  16-step  and  37-step  mechanisms  obtained  by  our  in- 
house  code  and  CHEMKIN  are  compared. 
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Temperature  (K) 


Figure  22:  Mole  Fraction  of  H2,  02,  OFI,  and  F120 
vs.  Time  for  the  16-step  Skeletal  Mechanism 


Figure  23:  Mole  Fraction  of  H2,  02,  OFI,  and  H20 
vs.  Time  for  the  37-step  Detailed  Mechanism 


Figure  20:  Temperature  vs.  Time  for  the  16-step 
Skeletal  Mechanism  and  Different  Values  of 
Equivalence  Ratio 


Figure  21:  Temperature  vs.  Time  for  the  37-step 
Detailed  Mechanism  and  Different  Values  of 
Equivalence  Ratio 
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Figure  24:  Mixture  Density  vs.  Time  for  the  16-step  Skeletal  and  the  37-step  Detailed  Mechanisms  (9  =  1.0) 

In  Figure  25,  the  Cp  of  the  mixture  has  been  depicted  versus  time  for  different  mechanisms. 


Figure  25:  Mixture  of  Cp  vs.  Time  for  Different  Mechanisms,  (9  =  1.0) 

As  mentioned  before,  there  are  three  different  methods  for  calculating  the  molecular  viscosity 
and  thermal  conductivity  of  the  species.  In  Figure  26  and  Figure  27,  these  three  methods  have 
been  compared  for  the  37-step  detailed  mechanism  and  equivalence  ratio  of  unity. 
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Figure  26:  The  Viscosity  of  H2,  02,  H20,  and  OH  Calculated  from  Different  Methods  for  the  37-step 

Detailed  Mechanism  vs.  Time,  (cp  =  1.0) 
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Figure  27:  The  Conductivity  of  H2,  02,  H20,  and  OH  Calculated  from  Different  Methods  for  the  37-step 

Detailed  Mechanism  vs.  Vime,  (cp  =  1.0) 


The  other  important  transport  property  is  the  binary  diffusion  coefficient  between  species  i  and  j. 
For  the  mechanisms  that  contain  9  species,  there  exist  45  distinct  binary  diffusion  coefficients. 
Figure  28  shows  the  trend  of  diffusion  coefficients  between  FF  and  other  species  for  the  37-step 
detailed  or  the  16-step  skeletal  mechanisms. 
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Figure  28:  The  Binary  Diffusion  Coefficients  between  H2  and  Other  Species  Multiplied  by  Pressure  vs. 

Temperature 

As  mentioned  before,  there  are  different  approaches  and  approximations  for  calculating  the 
mixture  viscosity  and  conductivity.  To  compare  these  different  approximations,  we  examine  the 
performance  of  all  methods  used  for  calculating  the  mixture  viscosity  and  conductivity  for  the 
same  problem,  i.e.  using  the  37-step  detailed  mechanism  for  the  equivalence  ratio  of  unity.  In 
this  problem  the  properties  of  the  individual  components  are  calculated  using  the  NASA 
formulas.  According  to  Figure  29,  the  third  approximation  yields  values  for  the  mixture 
viscosity  that  are  near  the  prediction  of  the  theoretical  method. 


Tme  fs'i 

Figure  29:  The  Mixture  Viscosity  Calculated  from  Different  Methods  vs.  Time  for  the  37-step  Detailed 

Mechanism,  (cp  =  1.0) 

As  shown  in  Figure  30,  the  first  approximate  method  has  the  best  prediction  for  the  mixture 
thermal  conductivity  compared  to  the  exact  method.  As  a  conclusion,  we  can  use  the  third  and 
first  approximate  methods  to  calculate  the  mixture  viscosity  and  thermal  conductivity, 
respectively,  instead  of  the  theoretical  methods  without  facing  any  significant  error.  It  must  be 
noted  that  if  we  use  this  combination  of  approximate  methods,  similar  to  what  has  been  done  in 
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Figure  29  and  Figure  30,  then  the  computational  cost  would  be  only  one-fifth  of  that  in  the  case 
in  which  the  exact  methods  are  used  for  the  calculation  of  mixture  properties. 


Time  (s) 


Figure  30:  The  Mixture  Thermal  Conductivity  Calculated  from  Different  methods  vs.  Time  for  the  37-step 

Detailed  Mechanism,  (cp  =  1.0) 


Finally,  it  is  useful  to  note  that,  if  one  uses  either  NASA  formula  for  the  transport  properties  of 
the  individual  components  or  the  third-order  polynomials  fitted  to  CHEMKIN  data,  one  obtains 
similar  results  for  the  mixture  properties.  This  is  confirmed  in  Figure  3 1  where  the  different 
methods  for  the  calculation  of  mixture  properties  have  been  used  for  the  hydrogen  combustion 
simulations  with  the  37-step  detailed  mechanism. 

For  the  simulations  of  other  (hydrocarbon)  fuels  such  as  ethylene  a  procedure  similar  to  that 
described  above  for  the  hydrogen  may  be  used.  However,  the  time  step  for  some  of  the  reaction 
steps  for  hydrocarbon  fuels  might  become  much  smaller  than  those  of  the  flow  and  mixing.  This 
requires  the  use  of  implicit  methods  and  ODE  solvers.  Additionally,  the  reaction  mechanisms 
for  hydrocarbons  have  to  be  reduced.  We  will  use  the  latest  reduced  mechanisms  for  these  fuels. 


Figure  31:  (a)  Mixture  Viscosity  and  (b)  Thermal  Conductivity  vs.  Time  for  the  37-step  Detailed  Mechanism 
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4.4  Combustion  of  Hydrogen  in  Spatially  Developing  Mixing  Layers 

The  spatially  developing,  reacting  subsonic  and  supersonic  mixing  layers  have  been  simulated 
with  the  hybrid  LES/FMDF  methodology  and  reaction  models  and  property  relations  described 
in  the  previous  section.  The  simulated  mixing  layer  is  composed  of  non-premixed  streams  of 
diluted  hydrogen  and  air.  The  chemical  reaction  of  hydrogen  and  air  is  modeled  with  the  37-step 
detailed  mechanism.  The  filtered  compressible  Navier-Stokes  equations  are  solved  with  high 
order  Eulerian  finite-difference  (FD)  methods.  The  hydrogen-air  mixing  and  combustion  are 
computed  by  solving  the  FMDF  equation  via  the  Lagrangian  Monte  Carlo  particle  method.  The 
consistency  between  the  LES-FD  and  FMDF-MC  methods  for  the  fuel  mass  fraction  and 
temperature  is  investigated  in  details  below. 

In  LES/FMDF  simulation  of  the  spatially  developing,  reacting  mixing  layer,  the  governing 
equations  are  solved  using  4th  order  spatial  compact  scheme  with  3ld  order  Runge-Kutta  method 
for  time  marching.  The  velocity  and  pressure  fields  are  calculated  using  the  finite  difference 
approach,  while  the  mass  fractions  of  all  nine  species  as  well  as  the  reaction  terms  are  computed 
using  the  Lagrangian  FMDF  solver.  Temperature  is  found  from  both  FD  and  FMDF  solvers  in 
order  to  investigate  the  consistency  between  the  predictions  of  the  Eulerian  and  Lagrangian 
approaches. 

In  the  simulated  subsonic  and  supersonic  mixing  layer  streams  are  composed  of  cold  diluted 
hydrogen  (2TF+N2)  in  one  stream  and  hot  air  (02+3.76N2)  in  another  stream.  After  initialization 
of  the  flow,  the  compressible  gas  dynamic  equations  are  solved  until  the  mixing  layer  is 
developed.  For  this  calculation,  the  distributions  of  fuel  and  oxidizer  are  also  obtained  at  each 
time  step  by  solving  a  filtered  passive  scalar  equation.  After  this  stage,  the  Monte  Carlo  particles 
are  distributed  in  the  domain  and  the  reacting  LES/FMDF  simulations  are  conducted.  At  the 
second  stage  all  mass  fractions  are  obtained  via  the  FMDF  method.  In  all  simulations,  six 
particles  per  grid  are  used.  The  transport  properties  such  as  the  viscosity  and  the  thermal 
conductivity  and  also  the  specific  heats  of  the  species  are  calculated  using  the  polynomials,  as 
explained  in  details  in  the  Subsection  3.4.  The  mixture  properties  are  calculated  using  relations 
(1-10),  (1-24),  and  (1-25).  In  all  simulated  cases,  molecular  Prandtl  and  Schmidt  numbers  are 
equal  to  0.72  while  the  turbulent  value  of  these  two  numbers  is  chosen  to  be  0.75.  In  addition,  the 
mixing  layer  thickness  at  inlet  is  2.5X10‘4  m. 

4.4.1  Subsonic  Case 

To  test  the  LES/FMDF  methodology  and  to  establish  the  reliability  and  the  accuracy  of  its 
numerical  solution  method,  we  have  first  simulated  subsonic  mixing  layers  with  and  without 
combustion.  Supersonic  mixing  layer  are  considered  next.  The  mixing  layer  inflow  and  initial 
conditions  are  shown  schematically  in  Figure  32.  The  tangent  hyperbolic  profiles  for  the 
streamwise  velocity  and  mass  fractions  of  fuel  and  oxidizer  are  used  to  initialize  the  domain.  The 
fuel  stream  has  a  lower  speed  and  temperature  than  the  oxidizer  stream. 
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Figure  32:  Schematic  View  of  Subsonic  Mixing  Layer  and  Initial  Conditions 

The  scatter  plots  of  fuel  mass  fraction  and  temperature,  obtained  from  LES-FD  and  FMDF-MC 
data,  for  the  nonreacting  subsonic  mixing  layer  are  shown  in  Figure  33.  There  is  a  good 
correlation  between  LES-FD  and  FMDF-MC  results,  indicating  the  accuracy  of  both  components 
of  the  hybrid  LES/FMDF  flow  solver  for  the  simulated  nonreacting  mixing  layer. 


U=300  m/s,  T=1100  K,  M=0.465 


U=100  m/s,  T=400  K,  M=0.15 


Figure  33:  Scatter  Plots  of  (a)  Hydrogen  Mass  Fraction,  and  (b)  Temperature,  Obtained  by  LES-FD  and 


FMDF-MC 


Figure  34  illustrates  the  snapshots  of  the  mole  fractions  of  the  nine  species  and  the  released  heat 
(J/Kg.s)  in  the  reacting  case.  In  this  simulation,  Smagorinsky  eddy  viscosity  coefficient  Cs  is 
0.028  and  the  filter  width  is  twice  as  the  computational  grid  size.  The  coefficient  C,P  in  the 
mixing  frequency  term  (Eq.  (62))  has  been  considered  to  be  3.5.  According  to  these  figures,  the 
mole  fraction  of  H02  and  H2C>2  are  small  compared  to  other  species.  But  contrary  to  the 
presumption  of  the  reduced  mechanisms,  mole  fraction  of  OH  is  comparable  to  those  of  H  and  O. 
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Figure  34:  Contours  of  Instantaneous  Species  Mole  Fractions  and  Heat  of  Reaction  (J/Kg.s) 

For  the  conditions  mentioned  above,  the  temperature  and  mass  fraction  of  hydrogen  predicted  by 
LES-FD  are  depicted  in  Figure  35  versus  the  corresponding  values  calculated  by  FMDF-MC. 
According  to  these  scatter  plots,  the  correlations  between  the  results  of  LES-FD  and  FMDF-MC 
are  very  high,  indicating  that  the  numerical  solution  of  LES/FMDF  equations  is  very  accurate. 
Figure  36  shows  the  distribution  of  the  temperature  calculated  by  LES-FD  and  FMDF-MC 
versus  the  mixture  fraction.  The  overall  trends  of  both  plots  as  well  as  the  upper  and  lower  limits 
are  the  same. 
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Figure  35:  Scatter  Plots  of  (a)  Hydrogen  Mass  Fraction,  and  (b)  Temperature,  Obtained  by  LES-FD  and 


FMDF-MC 


Figure  36:  Scatter  Plots  of  Temperature  Predicted 


(a)  LES-FD,  (b)  FMDF-MC  vs.  Mixture  Fraction 


In  order  to  understand  the  effect  of  mixing  frequency  on  the  consistency,  several  values  of 
coefficient  C,;,  were  examined  while  keeping  other  conditions  and  parameters  constant.  Figure  37 
and  Figure  38  show  the  scatter  plots  of  temperature  and  mass  fraction  of  FL  predicted  by  LES- 
FD  and  FMDF-MC  approaches  for  three  different  values  of  2.0,  3.5,  and,  5.0  of  Cv,  respectively. 
In  all  three  cases,  Smagorinsky  coefficient  is  0.028  and  the  filter  width  is  twice  as  the  FD  grid 
size.  According  to  these  graphs,  as  coefficient  C,,,  increases,  the  consistency  between  the  results 
is  slightly  improved.  Nevertheless,  the  consistency  between  LES-FD  and  FMDF-MC  results 
remains  to  be  very  high  in  all  cases.  This  again  indicates  the  reliability  of  the  hybrid  LES/FMDF 
methodology  and  its  numerical  solution.  It  should  be  noted  here  that  for  the  case  with 
combustion,  LES-FD  results  are  computed  by  using  the  reaction  source/sink  terms  obtained  from 
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FMDF  and  MC  particles.  This  is  only  possible  in  our  hybrid  LES/FMDF  solver  since  the 
reaction  is  closed  in  the  FMDF  formulation.  In  other  LES  models  the  very  nonlinear  and 
complex  SGS  reaction  terms  have  to  be  modeled! 


Figure  37:  Scatter  Plots  of  Hydrogen  Mass  Fraction  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  C(p=2.0  and 

(b)  C(p=3.5,  (c)  C(p=5.0  Conditions 


Figure  38:  Scatter  Plots  of  Temperature  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  C(p=2.0  and  (b) 

C(p=3.5,  (c)  C(p=5.0  Conditions 


Another  factor  that  might  affect  the  consistency  is  the  models  used  for  the  subgrid  scale  stress 
and  scalar  flux  models.  Ideally  we  do  not  expect  the  models  to  affect  the  consistency.  However, 
in  the  performed  simulations,  the  value  of  the  Smagorinsky  coefficient  and  the  filter  width  are 
significant.  The  increment  in  value  of  any  of  these  two  variables  leads  to  the  greater  value  of 
eddy  viscosity.  In  Figure  39(b)  and  Figure  40(b),  the  value  of  eddy  viscosity  is  four  times  greater 
than  that  of  case  (a).  For  both  simulations  Cv  is  3.5  and  other  conditions  are  the  same.  According 
to  these  graphs,  as  the  eddy  viscosity  increases,  the  consistency  is  slightly  improved. 
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Figure  39:  Scatter  Plots  of  Hydrogen  Mass  Fraction  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  A  =Grid 

Size  and  (b)  ^  =2xGrid  Size  Filters 


Figure  40:  Scatter  Plots  of  Temperature  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  ^  =Grid  Size  and  (b) 

A  =2  x  Grid  Size  Filters 


4.4.2  Supersonic  Case 

The  configuration  of  simulated  supersonic  mixing  layer  has  been  illustrated  in  Figure  41.  For 
this  set  of  conditions,  the  convective  velocity  and  Mach  number  are  equal  to  924  m/s  and  0.20, 
respectively.  Similarly  to  the  subsonic  case,  the  fuel  stream  has  a  lower  speed  and  temperature 
than  the  oxidizer  one. 
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U=1057  m/s.  T=1 100  K.  M=1.6 


Figure  41:  Schematic  View  of  Supersonic  Mixing  Layer  and  Initial  Conditions 


Figure  42  shows  the  scatter  plots  of  H2  mass  fraction  and  temperature  calculated  by  LES-FD  and 
FMDF-MC  methods  for  the  nonreacting  supersonic  mixing  layer.  In  this  simulation,  C(/)  is  equal 
to  3.5  and  the  Smagorinsky  coefficient  Cs  is  0.028. 
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Figure  42:  Scatter  Plots  of  (a)  Flydrogen  Mass  Fraction,  and  (b)  Temperature,  Predicted  by  LES-FD  and 

FMDF-MC 


The  snapshots  of  species  mole  fractions  as  well  as  the  combustion  heat  release  (J/Kg.s)  in  the 
reacting  case  are  presented  in  Figure  43.  In  this  simulation,  the  Smagorinsky  eddy  viscosity 
coefficient  is  0.028  and  the  filter  width  is  twice  as  the  computational  grid  size.  The  coefficient  Cf 
in  the  mixing  frequency  term  has  been  considered  to  be  3.5. 
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Figure  43:  Contours  of  Instantaneous  Species  Mole  Fractions  and  Heat  of  Reaction  (J/Kg.s) 


For  the  above-mentioned  conditions,  scatter  plots  of  temperature  and  mass  fraction  of  H2, 
obtained  by  LES-FD  and  FMDF-MC  methods  are  shown  in  Figure  44.  The  scattering  in  these 
data  is  less  than  those  shown  in  the  subsonic  flow.  Specifically,  the  mass  fractions  of  the  fuel 
calculated  by  Eulerian  and  Lagrangian  equations  are  highly  correlated  and  fully  consistent.  The 
instantaneous  values  of  the  temperature  obtained  by  LES-FD  and  FMDF-MC  in  the  whole 
domain  are  plotted  versus  mixture  fraction  in  Figure  45.  These  plots  again  show  the  very  good 
consistency  between  the  predictions  of  two  methods. 


Figure  44:  Scatter  Plots  of  (a)  Hydrogen  Mass  Fraction,  and  (b)  Temperature,  Predicted  by  LES-FD  and 

FMDF-MC 
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Figure  45:  Scatter  Plots  of  Temperature  Predicted  by  (a)  LES-FD  and  (b)  FMDF-MC  vs.  Mixture  Fraction 

In  order  to  understand  the  effect  of  Cf  on  the  consistency,  all  conditions  are  kept  the  same  but 
three  different  values  of  Cv  are  examined.  According  to  Figure  46  and  Figure  47,  the  consistency 
between  the  LES-FD  and  FMDF-MC  results  remain  unchanged  as  C,P  varies.  But  the  consistency 
between  the  temperatures  is  improved  slightly  as  Cf  is  increased. 


Figure  46:  Scatter  Plots  of  Hydrogen  Mass  Fraction  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  C(p=2.0,  (b) 

C(p=3.5  and  (c)  C(p=5.0  Conditions 


Figure  47:  Scatter  Plots  of  Temperature  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  C(p=2.0,  (b)  C(p=3.5 

and  (c)  C(p=5.0  Conditions 
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The  eddy  viscosity  has  a  significant  effect  on  the  consistency.  According  to  Figure  48  and  Figure 
49,  when  the  value  of  C(p  is  fixed  to  3.5,  as  the  filter  width  is  doubled,  the  fuel  mass  fraction  and 
temperature  become  more  consistent.  The  doubling  of  the  filter  width  can  be  viewed  as 
increasing  the  Smagorinsky  coefficient  Cs. 


Figure  48:  Scatter  Plots  of  Hydrogen  Mass  Fraction  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  ^  =Grid 

Size  and  (b)  ^  =2xGrid  Size  Filters 
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Figure  49:  Scatter  Plots  of  Hydrogen  Mass  Fraction  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  ^  =Grid 

Size  and  (b)  ^  =2xGrid  Size  Filters 

The  last  simulated  flow  considered  in  this  report  is  a  three-dimensional  spatially  developing 
supersonic  reacting  mixing  layer.  Due  to  the  large  number  of  computational  cells  and  MC 
particles,  a  parallel  processing  technique  for  the  hybrid  LES-FD/FMDF-MC  methodology  was 
used.  The  parallelization  method  is  similar  to  the  second  method  described  in  reference  [22].  The 
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initial  and  inflow  conditions  of  the  two  initially  non-premixed  fuel  and  oxidizer  streams  are 
shown  schematically  in  Figure  50.  For  this  set  of  conditions,  the  convective  velocity  and  Mach 
number  are  872.4  m/s  and  0.277,  respectively.  The  consistency  between  the  LES-FD  and  FMDF- 
MC  values  of  the  fuel  mass  fraction  and  temperature  in  the  simulated  3D  nonreacting  mixing 
layer  is  demonstrated  in  Figure  5 1 . 
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Figure  50:  Schematic  of  Supersonic  Mixing  Layer  and  Initial  Flow  Conditions 
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Figure  51:  Scatter  Plots  of  (a)  Hydrogen  Mass  Fraction,  and  (b)  Temperature  Predicted  by  LES-FD  and 


FMDF-MC  Methods 


Figure  52  shows  the  iso-surface  contours  of  temperature  predicted  by  LES-FD  and  FMDF-MC. 
The  high  temperature  zone  (red  colored)  between  the  fuel  and  oxidizer  streams  has  been  marked 
as  the  reaction  zone.  The  snapshots  of  species  mole  fractions  as  well  as  the  combustion  heat 
release  (J/Kg.s)  are  presented  in  Figure  53.  The  scatter  plots  of  temperature  as  predicted  by  LES- 
FD  and  FMDF-MC  versus  mixture  fraction  are  shown  in  Figure  54.  According  to  these  figures, 
the  overall  behavior  of  temperatures  obtained  by  LES-FD  and  FMDF-MC  are  the  same,  but  the 
results  of  LES-FD  are  more  scattered.  This  scattering  is  more  visible  along  the  line  on  the  right 
side  of  the  peak  temperature.  This  shows  the  weakness  of  the  FD  approach  when  dealing  with 
the  sharp  gradients.  The  overshoot  and  undershoot  in  the  LES-FD  results  are  due  to  limited  grid 
resolutions.  Contrary  to  LES-FD  results,  the  FMDF-MC  method  does  not  have  any  problem  in 
capturing  sharp  gradients  in  flow  variables,  thus  presents  a  smooth  temperature/mixture  fraction 
plots. 
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Figure  52:  Iso-Surface  Contours  of  Temperature  as  Predicted  by  (a)  LES-FD,  and  (b)  FMDF-MC  Methods 
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Figure  53:  Iso-surface  Contours  of  Instantaneous  Species  Mole  Fractions  and  Heat  of  Reaction  (J/Kg.s) 
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Figure  54:  Scatter  Plots  of  Temperature  Predicted  by  (a)  LES-FD  and  (b)  FMDF-MC  vs.  Mixture  Fraction 


To  study  the  effect  of  mixing  term  on  the  consistency,  two  values  of  2.5  and  3.5  were  used  for 
the  coefficient  Cf  while  other  conditions  were  kept  the  same.  Figure  55  and  Figure  56  show  the 
scatter  plots  of  mass  fraction  of  H2  and  temperature  for  these  two  cases.  Similar  to  the  planar 
reacting  subsonic  and  supersonic  mixing  layers,  as  the  value  of  Cf  is  increased,  the  FD  and  MC 
predictions  become  slightly  more  correlated.  But  there  is  an  optimum  value  for  C(p  that  was 
found  to  be  around  3.5  for  which  the  highest  consistency  is  achieved. 
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Figure  56:  Scatter  Plots  of  Temperature  Predicted  by  LES-FD  and  FMDF-MC  for  (a)  C(p=2.5  and  (b)  Cq>=3.5 

Conditions 

4.5  Implementation  of  FMDF  into  DSEM 

Several  simulations  were  perfonned  on  the  free  shear  flow  case  to  validate  the  particle-tracking 
algorithm  as  well  as  the  calculation  and  comparison  of  scalar  density.  In  the  first  test  case  the 
flow  is  periodically  forced  in  the  vertical  direction  at  the  inlet  to  accelerate  the  formation  of 
vortical  structures.  At  each  time  step  particles  are  seeded  at  predetermined  locations  along  the 
inlet  and  are  assigned  the  same  velocity  as  the  local  fluid  velocity.  On  the  contrary,  the  second 
flow  does  not  include  forcing  and  particles  are  seeded  at  a  constant  mass  flow  rate  in  a  uniformly 
random  fashion  across  the  inlet.  The  second  case  allows  us  to  test  the  consistency  between  the 
Eulerian  and  Lagrangian  approximations. 

4.5.1  Effect  of  Inertia  on  Fluid  Particles 


500  1000  1500  2000 

f 


(a) 

_Correlation 
.  Coefficient=  0.9978 


500  1000  1500  2000 

f 


(b) 

.Correlation 
Coefficient=  0.9982 


As  a  first  test,  we  considered  both  massless  and  inertial  particles  to  assess  the  prediction  of  the 
code  against  the  well-known  physical  behaviors.  For  inertial  particles  we  used  the  full  particle 
equations  of  motion,  which  includes  the  Lagrangian  equations  for  both  position  and  velocity  of 
the  particle. 

Three  cases  were  simulated  with  varying  Stokes  numbers  for  the  particles.  The  Stokes  number  is 
defined  as  the  ratio  of  the  particle  time  constant  to  the  carrier  phase  time  constant 


St  = 


(105) 


Here  rp  is  the  characteristic  time  attributed  to  the  particle,  otherwise  referred  to  as  response  time, 
Uf  is  the  characteristic  velocity  of  the  flow,  typically  taken  as  the  inlet  velocity  in  a  shear  flow, 
and  L*f  is  the  characteristic  length.  The  particle  routines  in  the  DSEM  code  accept  rp  as  a 
parameter,  and  the  fluid  characteristics  are  based  on  boundary/initial  conditions.  Note  that  for  a 
massless  fluid  particle  St=0. 
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For  the  case  of  a  high  Stokes  number,  St=2.73,  the  inertial  particles  characteristic  time  is  greater 
than  the  largest  time  scales  of  the  carrier  phase.  Therefore,  the  particles  do  not  respond  to  small 
scales  of  the  flow  and  escape  the  vortical  structures  formed  by  the  shear  layer.  Such  an  effect  is 
demonstrated  in  Figure  57,  where  the  initial  forcing  condition  in  the  y-direction  applied  at  the 
inlet  causes  the  particles  to  spread  in  the  y-direction  of  the  domain.  Note  that  the  figure  shows 
the  locations  of  all  the  particles  present  in  the  computational  domain  at  a  given  time. 
Alternatively,  for  St=1.0  the  characteristic  times  of  the  fluid  and  particles  are  comparable. 
Although  the  particles  carry  inertia,  the  flow  field  affects  their  behavior.  Such  an  effect  is  seen  in 
Figure  58,  where  the  profile  of  shear  flow  may  be  seen  in  the  particles  locations.  Finally,  Figure 
59  shows  that  massless  fluid  particles  follow  vortical  motions  of  the  flow,  and  their  locations 
resemble  the  Eulerian  flow  field.  Such  particles  are  necessary  for  the  FMDF  method,  and  it  may 
be  seen  in  Figure  59  that  they  behave  as  expected. 
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Figure  58:  Locations  of  all  the  Particles  with  St  =  1  at  a  Given  Time 


Figure  59:  Locations  of  all  the  Massless  (Fluid)  Particles  (St=0)  at  a  Given  Time 
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4.5.2  Comparison  of  Eulerian  and  Lagrangian  Approximations 


A  subsonic  free  shear  flow  was  selected  as  the  test  case  for  the  hybridized  DSEM/FMDF  code. 
The  DSEM  and  ensemble  grids  may  be  seen  in  Figure  3  and  Figure  4,  respectively.  Initial 
conditions  for  velocity  and  density  were  imposed  on  the  carrier  phase  as  follows: 


u  = 


t/i  +  U2 
*  -  + 


tanh 


v  —  0 


P  - 


P1+P2 

2 


+ 


Pi  ~  P 2 


tanh 


where  U1rdnd  U2  are  the  fluid  inlet  velocities  in  the  top  and  bottom  layers,  respectively,  and  Sw  is 
the  vorticity  thickness.  The  above  conditions  were  then  interpolated  to  the  particle  locations  as 
seen  in  Figure  60  for  density.  Similarly,  the  inlet  conditions  were  tangentially  hyperbolic  along 
the  y-axis  with  respect  to  velocity  and  density.  Unlike  the  previous  section,  no  velocity  or 
forcing  was  imposed  along  the  y-direction.  The  Lagrangian  boundary  conditions  were  imposed 
by  taking  advantage  of  the  matching  condition  described  in  the  section  above.  The  set  values  for 
the  constants  were 
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Figure  60:  Shear  Flow  Initial  Condition:  Lagrangian  Particles  and  Eulerian  Density  Contours 
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Figure  61:  Fluid  Density  Profile  at  x=0.30,  Calculated  by  Eulerian  and  Lagrangian  Methods 


Figure  61  shows  consistency  between  our  Eulerian  and  Lagrangian  methods.  In  this  figure,  we 
depict  the  cross-stream  density  sampled  at  a  distance  of  x=0.30  along  the  stream.  It  is  evident 
that  the  Lagrangian  (FMDL)  prediction  is  consistent  to  the  Eulerian  (DSEM)  prediction. 

The  oscillations  seen  in  Ligure  61  may  be  attributed  to  statistical  error.  Lrom  a  statistical 
standpoint  it  would  be  ideal  to  have  infinitely  many  particles  in  an  infinitely  small  domain. 
However  due  to  computational  overhead  it  is  necessary  to  select  a  sufficient  number  of  particles 
with  sufficiently  small  domains  as  to  ensure  statistical  accuracy  while  minimizing  dispersive 
error.  Such  parameters  cannot  be  predicted  a  priori ;  however  it  is  typically  considered  that  for 
such  a  simulation  40  particles  per  ensemble  domain  are  sufficient.  Yet  in  areas  of  high 
variability,  the  impact  of  a  single  particle  entering  or  leaving  a  domain  may  manifest  itself  as  an 
oscillation  in  the  density.  In  order  to  minimize  such  oscillation  it  is  necessary  to  adjust  the 
weighting  of  said  particles  so  that  an  individual  particle  carries  less  contribution  to  the  overall 
calculation  of  density  within  the  ensemble  domain  and  a  larger  sample  must  be  taken  [17].  To 
employ  such  a  method,  it  would  be  necessary  to  first  parallelize  the  DSEM  code  so  that  more 
particles  may  be  computed  in  a  timely  manner.  This  will  be  done  in  Phase  II. 
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5.0  Conclusions 


The  ultimate  goal  of  this  project  is  to  develop  and  validate  a  novel  numerical  code  for  simulation 
of  high  speed,  reacting  flows  based  on  several  innovative  methods.  The  work  during  Phase  I  was 
devoted  to  providing  proof-of-concept  for  each  of  the  following  three  components  of  the  study: 

1 .  Implementation  of  the  EV  method  in  the  DSEM  code 

2.  Development  of  FMDF  for  high-speed  flows 

3.  Implementation  of  FMDF  in  DSEM 

In  Phase  I,  in-house  ID  and  2D  spectral  element  codes  were  developed  to  investigate  the  idea  of 
using  an  artificial  dissipative  term  in  conservation  laws  to  avoid  numerical  oscillations  due  to 
high-order  approximations  in  conjunction  with  abrupt  changes  in  the  solution  (physical  shocks). 
Increasing  the  magnitude  of  the  local  viscosity  proportional  to  the  rate  of  local  entropy 
production  should  not  affect  the  solution  where  the  flow  is  smooth  while  efficiently  removing 
spurious  oscillations. 

In  ID  simulations,  explicit  filtering  was  also  tested  as  an  alternative  remedy  and  the  results  of  the 
two  approaches  were  compared  for  two  model  problems.  For  the  linear  advection  problem,  the 
results  showed  that,  filtering  is  effective  in  removing  high-frequency  oscillations  but  does  not 
produce  a  completely  oscillation-free  solution.  However,  this  was  not  the  case  for  the  inviscid 
Burgers  problem.  The  ID  results  obtained  by  entropy  viscosity  were  promising  and  even  for  very 
sharp  discontinuities  this  method  was  able  to  produce  relatively  smooth  solutions  and  remove 
major  instabilities  in  the  solution. 

In  2D  simulations,  where  we  focused  on  Euler  equations,  a  special  attention  was  paid  to 
boundary  and  interface  treatment  in  the  2D,  DSEM-EV  code.  This  effort  led  to  successful 
implementation  of  the  EV  method  in  2D.  Specifically,  the  interface  issue  that  was  observed  in 
our  earlier  results  was  completely  addressed  and  resolved. 

Several  test  cases  were  simulated  with  the  DSEM-EV  code  and  an  effort  was  made  to  quantify 
the  magnitude  of  the  error  in  each  case  against  a  (high-resolution)  reference  case.  It  was  shown 
that,  the  error  decreases  by  increasing  either  the  number  of  elements  or  the  polynomial  order.  For 
the  same  total  number  of  nodes,  using  higher  polynomial  order  decreases  the  magnitude  of  the 
error  more  effectively.  Consistency  and  stability  of  the  code  was  demonstrated  and  the  code  is 
ready  to  be  used  for  more  complex  geometries  in  2D. 

Development  of  FMDF  for  high-speed  flows  started  with  developing  an  in-house  reaction 
module  for  numerical  simulation  of  hydrogen  in  air  combustion  with  different  chemical  kinetics 
mechanisms.  Different  methods  for  calculating  the  species  and  mixture  properties  were 
investigated,  and  efficient  ways  to  calculate  the  thermodynamic  and  transport  properties  were 
chosen.  The  results  obtained  by  our  reaction  module  for  a  premixed  stirred  reactor  were 
compared  with  those  of  CHEMKIN. 

The  EES-FMDF  method  was  used  together  with  the  37-step  detailed  kinetics  mechanism  for  2D 
and  3D  simulations  of  hydrogen-air  combustion  in  planar  subsonic  and  supersonic  mixing  layers. 
The  consistency  between  the  species  mass  fractions  and  temperatures  predicted  by  EES-FD  and 
FMDF-MC  methods  was  established.  The  effects  of  mixing  frequency  and  eddy  viscosity  on  the 
computed  results  and  the  consistency  were  reported.  The  EES-FMDF  method  is  ready  to  be  used 
for  more  complex  high  speed,  reacting  flows  involving  complex  chemical  kinetics. 
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Preliminary  work  started  in  Phase  I  for  implementation  of  FMDF  in  DSEM.  A  particle-tracking 
algorithm  suitable  for  use  with  FMDF  was  developed  and  tested  in  the  2D  DSEM  code. 
Furthermore,  with  the  addition  of  the  particle  density  and  weighting  routines,  we  were  able  to 
approximate  a  scalar,  fluid  density,  in  the  Lagrangian  frame  and  perform  a  consistency  check 
against  the  Eulerian  approximation.  We  showed  the  accuracy  of  our  FMDF  algorithm  to 
produce  suitable  results.  With  additional  work  and  the  implementation  of  additional  SDE’s,  it 
will  be  possible  to  solve  for  other  scalars  such  as  temperature  and  reactions. 
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AB 

CFD 

DPT 

DSEM 

EV 

EES 

FD 

FMDF 

gg 

gH 

LES 

lg 

MC 

PDE 

PDF 

SDE 

SGS 

St 


List  of  Acronyms,  Abbreviations,  and  Symbols 

Adams-Bashforth 
Computational  Fluid  Dynamics 
Discrete  Polynomial  Transform 
Discrete  Spectral  Element  Method 
Entropy  Viscosity 
Enabling  Energy  Systems 
Finite  Difference 
Filtered  Mass  Density  Function 
Gauss-Gauss 
Gauss-Lobatto-Lobatto 
Large-Eddy  Simulation 
Lobatto-Gauss 
Monte  Carlo 

Partial  Differential  Equation 
Probability  Density  Function 
Stochastic  Differential  Equation 
Sub-Grid  Scale 
Stokes  Number 
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